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1.  INTRODUCTION 

1.1.  Scope  of  Thesis 

Multi-level,  or  multi-grid,  methods  have  existed  for  more  than 
fifteen  years,  but  have  only  recently  become  popular.  Thus,  there 
remain  many  unanswered  theoretical  questions.  Several  of  these  are 
addressed  in  this  thesis,  the  major  one  being  the  computational  complexity 
of  multi-level  methods  for  locally  refined  finite  element  grids.  These  are 
grids  for  which  the  ratio  of  diameters  of  the  largest  and  smallest 
elements  is  very  large,  becoming  unbounded  on  the  family  of  discretizations 
of  interest.  The  principal  result  of  this  thesis  is  that  multi-level 
methods  can  approximately  solve  the  linear  algebraic  systems  for  such 
finite  element  discretizations  in  0(N)  operations  for  a  grid  with  N 
elements.  Such  asymptotically  optimal  complexity  bounds  have  been  given 
previously  for  multi-level  methods  on  non-locally-ref ined  grids,  but  not 
for  locally  refined  grids.  This  result  provides  theoretical  justification 
for  the  use  of  multi-level  methods  on  locally  refined  grids,  which  is 
becoming  increasingly  popular.  The  proof  of  this  result,  which  occupies 
almost  half  of  this  thesis,  is  also  of  interest  for  several  reasons. 

Perhaps  the  most  significant  of  these  is  that  the  rate  of  convergence 
shown  for  the  multi-level  iteration  used  depends  only  on  local  properties 
of  the  finite  element  spaces  involved.  Global  properties,  including 
elliptic  regularity,  are  never  invoked.  As  a  consequence  this  result 
provides  an  independent  verification  of  the  0(N)  complexity  of  multi-level 
methods  for  non-convex  domains,  established  previously  by  Bank  and  Dupont 
(1978).  It  also  opens  the  door  to  easily  computable  bounds  on  the 
convergence  rate  of  multi-level  methods  on  irregular  finite  element  grids. 
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In  all,  three  multi-level  algorithms  are  considered  in  detail 
in  this  thesis.  The  first,  analyzed  in  chapter  three,  is  similar  to 
methods  considered  by  Nicolaides  (1977)  and  by  Bank  and  Dupont  (1978), 
but  is  simpler  than  these  others.  It  is  shown  to  converge  in  0(N) 
operations  for  N-point  quasi-uniform  finite  element  grids  on  convex 
domains,  through  a  convergence  proof  similar  to  those  given  by  the  above 
authors.  This  proof  is  given  here  since  it  is  less  complex  than  an 
earlier  one  based  on  Fourier  analysis  of  the  truncation  error.  The 
algorithm  for  quasi-uniform  grids  is  treated  here  partly  as  background 
for  the  analysis  of  algorithms  for  locally  refined  grids  in  chapter  four. 
It  is  also  of  independent  interest  that  an  0(N)  bound  can  be  shown  for 
this  algorithm  even  though  it  is  simpler  than  those  in  the  literature. 

Two  algorithms  are  considered  in  chapter  four,  both  applicable 
to  locally  refined  grids  on  arbitrary  polygonal  domains.  The  first  of 
these  is  quite  similar  to  the  method  used  in  the  adaptive  finite  element 
code  of  Bank  and  Sherman  (1978).  It  is  the  closest  algorithm  to  theirs 
that  could  be  analyzed  by  the  theoretical  techniques  here.  This  method 
is  shown  to  converge  in  0(N)  operations  subject  to  a  condition  on  the 
rate  of  growth  of  the  dimensions  of  the  finite  element  spaces  used.  The 
necessary  growth  condition  requires  that  the  finer  finite  element  spaces 
used  by  the  multi-level  algorithm  have  far  more  elements  than  the 
coarser  spaces  used.  Such  a  condition  is  also  imposed  in  the  Bank  and 
Sherman  code  though  it  is  less  stringent  there  due  to  their  technique 
of  "level  compression,"  Bank  and  Sherman  (1978). 

Since  these  growth  conditions  can  be  quite  restrictive,  a 
second  algorithm  is  also  considered  in  chapter  four.  It  is  a  modification 
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of  the  first,  and  it  will  be  shown  to  converge  in  0(N)  operations  under 
a  considerably  weaker  growth  condition  on  the  dimensions  of  the  finite 
element  spaces  used.  Though  the  theoretical  interest  in  this 
modification  is  apparent,  the  author  is  at.  present  unsure  of  the 
practical  utility  of  this  second  algorithm.  One  pays  a  large  penalty 
for  the  weaker  growth  conditions  through  an  increase  in  the  hidden 
constant  in  the  0(N)  work  bound. 

1.2.  History  of  Multi-Level  Methods 

Multi-level  or  multi-grid  methods  are  usually  thought  of  as 
originating  in  the  work  of  Fedorenko  (1961,  1964),  although  related  ideas 
had  been  put  forth  earlier.  The  idea  of  using  several  discretizations 
in  the  solution  of  a  PDE  (partial  differential  equation)  is  actually 
quite  old.  Southwell  (1940)  suggested  using  the  solution  of  an  elliptic 
problem  on  a  coarse  finite  difference  grid  to  obtain  a  starting  value 
for  relaxation  on  a  finer  grid.  He  also  recommended  a  technique  called 
"block  relaxation"  which  is  closely  related  to  multi-level  methods. 

In  this  technique  the  person  doing  the  relaxation  calculation  (in  the 
era  before  computers)  would  find  a  region  of  the  finite  difference  grid 
where  the  residuals  were  all  of  the  same  sign.  A  constant  would  then 
be  added  to  the  value  of  the  trial  solution  throughout  this  region  to 
reduce  the  average  residual  on  this  region  as  much  as  possible.  As  a 
consequence,  after  this  correction  the  residuals  would  oscillate  in  sign 
on  this  region.  Southwell  observed  that  as  long  as  there  were  no  large 
areas  where  the  residual  did  not  alternate  in  sign,  relaxation  converged 
rapidly.  Of  course,  though  quite  effective,  such  a  heuristic  approach 
was  not  of  much  help  in  the  early  days  of  computers,  so  simpler  methods 


like  SOR  took  over. 
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The  idea  of  making  block  corrections  essentially  died  out  until 

Fedorenko  (1961)  observed  that  rather  than  making  corrections  in  isolated 

blocks  as  suggested  by  Southwell,  one  could  make  a  global  correction  by 

solving  a  so-called  coarse  grid  residual  problem.  For  example,  suppose 

that  in  iteratively  solving  Poisson's  equation  on  a  grid  with  mesh  size 

h  one  has  at  the  n-th  iteration  a  trial  vector  u^  £  G,  and  residual 

— n  h 

_r^  £  where  G^  is  the  vector  space  of  mesh  functions  on  the  grid. 

If  A^  is  the  discrete  Laplacian  being  used. 


h  h  .  h 
r  =  f  -  A,  u 
— n  —  h  — n 


where  _f  is  the  data  vector.  Now  let  G2h  be  the  vector  space  of  mesh 
functions  on  a  coextensive  coarser  grid  with  mesh  size  2h  and  suppose  we 
have  an  interpolation  mapping 


X2hS  G2h  *  Gh  * 

2h 

One  can  ask  for  a  correction  vector  v  £  G„,  such  that  the  improved 

—  zh 

solution 


h  h  _h  2h 

u  ...  =  u  +  I„,  v 
— n+1  — n  2h  — 

would  have  smaller  residual.  This  idea  has  been  considered  independently 

2h 

bv  others,  Wachspress  (1974).  It  is  attractive  since  v  lies  in  the 
lower  dimensional  space  G^  and  so  should  be  relatively  easy  to  compute. 

21i 

Fedorenko's  fundamental  contribution  was  to  note  that  v  could  be 
chosen  as 


(1.2.1) 


^2h  — 


2h 


x2h  rh 
h  — n 


,2h 


where  I,  :  G,  -*•  G„,  is  an  interpolation  mapping  and  A„,  is  a  discrete 
n  h  Zn  zn 

Poisson  operator  on  62^-  Since  (1.2.1)  is  just  a  discretization  of 

2  h  h 

Poisson's  equation  with  data  I,  r  , he  realized  that  (1.2.1)  could  be 

h  — n 
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approximately  solved  by  relaxation  and  the  recursive  use  of  corrections  on 
a  still  coarser  grid  G^. 

Fedorenko’s  theoretical  work  showed  that  for  a  model  problem 
multi-level  iteration  reduced  the  convergence  error  by  an  amount  per 
iteration  independent  of  the  mesh  size.  Thus, reducing  the  convergence 
error  to  a  magnitude  comparable  to  the  truncation  error  required 
0(N  log  N)  operations  on  an  N  point  grid,  the  same  as  ADI.  Shortly 
thereafter  Bakhvalov  (1966)  considered  beginning  the  solution  process 
on  a  very  coarse  grid,  using  the  solution  there  as  a  starting  value  for 
multi-level  iteration  on  a  slightly  finer  grid,  and  so  on,  gradually 
descending  down  to  the  level  of  refinement  desired.  This  idea, which 
was  probably  old  when  Southwell  was  young,  changes  the  work  bound  from 
0(N  log  N)  to  0(N).  No  other  current  method,  either  direct  or  iterative, 
achieves  this  asymptotically  optimal  work  bound. 

Besides  improving  the  asymptotic  work  bound  on  multi-level 
iteration,  Bakhvalov  extended  the  method  to  a  large  class  of  PDEs  on 
a  square  finite  difference  grid.  After  his  work,  almost  no  work  seems 
to  have  been  done  on  multi-level  methods  until  this  decade.  Then  several 
people  became  interested  in  this  method  including  Brandt,  Hackbusch  and 
Nicolaides.  Brandt  in  particular  is  largely  responsible  for  popularizing 
this  method  and  developing  it  into  a  practical  computational  tool.  The 
other  two  people  mentioned  have  been  working  more  on  the  theoretical 
questions  involved  in  this  method  and  on  generalization  to  finite  element 
grids. 

In  the  last  few  years  there  has  been  a  rapid  growth  of  interest 
in  multi-level  methods,  and  there  is  considerable  ongoing  research.  This 
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research  is  directed  both  at  the  theoretical  problems.  Bank  and  Dupont 
(1978),  Hackbusch  (1978),  Nicolaides  (1977),  and  at  developing  algorithms 
for  the  many  applications  areas.  Representative  work  of  the  latter 
type  may  be  found  in  Bank  and  Sherman  (1978),  Brandt  (1977),  Brandt, 

Dendy  and  Ruppel  (1978),  Nicolaides  (1975),  and  Poling  (1977).  If 
the  author  is  not  mistaken,  we  should  see  continued  research  in  this 
area  for  many  years  and  may  eventually  see  multi-level  methods  as 
the  method  of  choice  for  most  parabolic  and  elliptic  problems. 


1.3.  The  Finite  Element  Approach 

Design  and  analysis  of  multi-level  methods  can  be  done  from 
either  a  finite  difference  or  finite  element  point  of  view.  In  one 
sense  these  two  approaches  are  almost  the  same;  finite  element  methods 
are  a  special  class  of  finite  difference  method.  Nevertheless,  the  way 
people  think  about  finite  element  methods  and  the  kinds  of  grids  they 
use  for  these  methods  are  so  different  from  the  usual  way  people 
understand  and  use  finite  difference  methods  that  we  may  think  of  these 
classes  as  completely  separate.  This  thesis  is  devoted  to  analysis 
of  multi-level  methods  from  the  finite  element  point  of  view.  Several 
of  the  many  reasons  the  author  had  for  adopting  this  point  of  view  will 


be  described  in  this  section. 
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While  the  multi-level  method  originated  as  a  means  of  solving 
finite  difference  equations  and  was,  until  recently,  solidly  wedded  to 
the  finite  difference  approach,  much  of  the  current  work  on  multi-level 
methods  is  being  done  in  the  finite  element  context.  Starting  with 
Hackbusch  (1977)  and  Nlcolaides  (1977),  who  seem  to  have  been  the  first 
to  consider  multi-level  methods  as  a  means  of  solving  finite  element 
equations,  there  has  been  increasing  interest  in  this  approach.  The 
reasons  for  this  are  not  complicated.  Finite  element  theory  is  in  many 
respects  superior  to  finite  difference  theory.  A  great  many  problems 
in  the  finite  difference  approach  either  do  not  arise  in  the  finite 
element  approach  or  are  easily  handled  there.  Questions  of  stability, 
the  use  of  irregular  grids,  treatment  of  corner  singularities,  treatment 
of  curved  boundaries  and  construction  of  error  estimates  are  all  much 
simpler  in  the  finite  element  context. 

Aside  from  the  obvious  (to  the  author)  superiority  of  the 
finite  element  method,  there  are  special  reasons  unique  to  multi-level 
algorithms  for  doing  their  theory  in  the  finite  element  setting.  The 
first  of  these  is  that  the  finite  element  error  estimates  are  based  on 
approximation  results  rather  than  on  asymptotic  expansions.  What  this 
means  for  multi-level  methods  is  that,  unlike  some  of  the  early  finite 
difference  arguments  where  rapid  convergence  of  the  iteration  could  be 
shown  only  for  sufficiently  small  size,  no  such  proviso  is  ever 
needed  in  the  finite  element  context.  A  second,  more  important  reason 
is  that  multi-level  algorithms  always  involve  more  than  one 
discretization.  Thus  one  requires  mappings  from  the  space  of  grid 


functions  or  finite  element  functions  of  each  discretization  to  those  of 
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other  discretizations.  In  the  finite  difference  context  such  mappings 
are  done  by  interpolations  which  require  separate  error  analyses.  In 
the  finite  element  context  the  grids  can  be  designed  so  that  the 
finite  element  space  on  a  coarser  grid  will  be  a  subspace  of  the  finite 
element  spaces  on  finer  grids.  Then  the  mapping  from  a  coarser  grid  to 
a  finer  grid  will  be  just  the  natural  injection.  To  go  the  other  way, 
some  kind  of  projection  is  needed,  but  many  natural  ones  are  available. 
The  fact  that  different  finite  element  discretizations  can  be  related 
in  this  way  yields  a  great  simplification  of  the  theory  of  multi-level 
methods.  Unfortunately  only  the  theory  is  simplified,  not  the 
programming.  Rather  than  an  interpolation  formula,  one  has  a  similar 
formula  relating  the  basis  functions  of  one  finite  element  space  to 
those  of  another.  It  is  interesting  that  in  many  cases  these  change 
of  basis  formulas  arising  in  the  finite  element  setting  are  exactly 
the  same  as  the  interpolation  formulas  already  in  use  in  the  finite 
difference  setting. 

There  is  another  reason  why  studying  multi-level  convergence 
is  important  for  this  thesis.  In  chapter  four  we  treat  multi-level 
methods  for  locally  refined  grids.  While  it  may  be  possible  to  analyze 
multi-level  methods  for  locally  refined  finite  difference  grids  in 
special  cases,  it  seems  unlikely  that  an  analysis  as  general  as  that 
given  here  would  be  possible  outside  finite  element  theory.  For 
locally  refined  grids,  finite  element  theory  also  helps  in  coding. 

One  never  needs  to  choose  special  formulas  for  points  where  the  mesh 
size  changes  abruptly  or  worry  about  truncation  error.  There  are 
also  beginning  to  be  reliable  a  posteriori  error  estimates  for  finite 


elements  enabling  one  to  design  robust  adaptive  algorithms,  Babushka 
and  Rheinboldt  (1978).  In  short,  multi-level  methods,  adaptive 
refinement,  and  finite  element  theory  mesh  together  quite  nicely. 

The  main  limitation  of  the  finite  element  approach  is  that  it 
is  far  more  natural  and  powerful  for  self-adjoint  problems  than  for  the 
general  case.  Neither  finite  difference  nor  finite  element  methods  work 
very  well  for  non-self-adjoint  problems,  but  it  is  frequently  easier 
to  see  what  to  do  on  regular  finite  difference  grids  than  on  general 
finite  element  grids.  While  it  would  have  been  possible  to  extend 
many  of  the  results  of  this  thesis  to  the  non-self-adjoint  case  by 
adding  various  restrictions,  that  was  a  more  ambitious  project  than  the 
author  was  prepared  to  attempt.  If  one  takes  the  results  here  as 
indicative  of  the  kind  of  results  that  could  be  proved  in  the  general 
non-self-adjoint  case,  one  should  not  go  far  wrong. 
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2.  PRELIMINARIES 

In  this  chapter  we  collect  notations  and  basic  results  from 
finite  element  theory.  Since  this  material  is  completely  standard,  the 
reader  familiar  with  finite  elements  can  probably  skip  this  chapter 
entirely,  referring  back  to  it  later  as  required.  The  reader  less 
familiar  with  the  finite  element  approach  to  partial  differential 
'  equations  may  find  this  chapter  helpful. 

2.1.  Elliptic  Equations 

Throughout  this  thesis  we  will  be  concerned  with  solving 
self-adjoint  second  order  boundary-value  problems.  The  problem  of  most 
interest  will  be  the  Dirichlet  problem 
(2. 1.1a)  Lu  =  f  on  !! 

(2.1.1b)  u  =  0  on  , 

2 

where  0,  is  a  bounded  open  domain  in  the  plane  IR  with  piecewise  smooth 
boundary  9£2.  L  may  be  any  self-adjoint  second  order  operator,  for 
example  the  Laplacian,  or  more  generally  any  operator  in  divergence  form 

Lu = -V  •  (a  V  u)  +  bu 

_°o 

where  a,  b  are  C  functions  on  fi.  In  order  for  this  operator  to  be 
elliptic  on  fi,  we  require  that  there  be  constants  a,  a,  b  such  that 

a  1.  a  £  a 

0  <  b  <  b  j 

hold  throughout  ft.  We  will  also  be  considering  the  Neumann  problem 

(2.1.2a)  Lu  *  f  on  SI  ] 

(2.1.2b)  u  -  0  on  30  , 

n 

where  u^  is  che  derivative  of  u  in  the  direction  normal  to  the  boundary. 

i 


4 


« 


* 
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In  this  case  ellipticity  requires  that  b  be  positive  on  some  neighborhood 
in  fi.  Bank  and  Dupont  consider  also  the  singular  Neumann  problem,  b  =  0, 
but  we  exclude  it  here. 

The  finite  element  approach  to  boundary-value  problems  is  based 

in  part  on  the  functional  analytic  theory  of  partial  differential 

equations,  in  which  the  boundedness  or  invertibility  of  a  differential 

operator  is  expressed  in  terms  of  Sobolev  norms.  These  norms  are  a 

2 

natural  generalization  of  the  L  norm,  defined  for  non-negative  integers 
l  by 


ML  -  <  e  ll^ull2,)1'2 

0<  a  <£  L 


where  a  =  (a^,  a  )  is  a  multi- index  with  a^,  >  0  and 


l«l  ■  +  a2 


Da  =  3 


al  a2 
3x  9y 

This  definition  makes  sense  for  all  functions  on  ft  whose  distribution 

2  OO 

derivatives  up  to  order  l  exist  and  lie  in  L  .  The  completion  of  C  (Q) 
in  || -||  is  denoted  and  is  a  Hilbert  space  relative  to  the  inner 


product 


(u,v)^  =  /  (  £  (Dau)(Dav))  . 


0<Ja|<fc 

If  1  is  a  negative  integer  ||  •  ||  is  defined  by  duality. 


sup 

veH£(ft) 


v) 


C  • 


Here  and  throughout,  (• ,  •)  is  the  L  inner  product  while  |j  •  ||  with  no 

2 

subscript  denotes  the  L  norm.  Sobolev  spaces  can  also  be  defined  for 
arbitrary  real  l  either  through  the  interpolation  theory  of  Hilbert 


« 

| 


* 

| 
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spaces  or  by  Fourier  transforms.  Details  may  be  found  in  Oden  and  Reddy 
(1976). 

Equations  (2.1.1)  and  (2.1.2)  implicitly  require  the  solution 
u  to  be  twice  differentiable,  which  is  overly  restrictive  for  many  purposes. 
Consider  the  problem 

(2.1.3a)  Lu  =  f  on  Q 

(2.1.3b)  u  =  0  on  T 

(2.1. 3c)  u  =  0  on  T  , 

n  2 

where  u  which  includes  both  (2.1.1)  and  (2.1.2)  as  special 

1  1 

cases.  Let  H  (fi)  be  the  subspace  of  H  (0)  satisfying  the  (essential) 

E 

Dirichlet  boundary  condition  (2.1.3b)  but  otherwise  unrestricted.  If 

(2.1.3)  is  satisfied  then 

(Lu,v)  =  (f,v),  ve  H^(Q)  . 

The  notation  here  means  that  equality  holds  for  all  test  functions  v  in 
H^(f2).  Defining  the  bilinear  form 

a(u,v)  =  (Lu,v)  , 

we  can  ask  for  the  element  u  in  H^(fi)  such  that 

£ 

(2.1.4)  a(u,v)  =  (f.v),  ve  hJ(£1)  . 

It  is  only  necessary  that  u  lie  in  H^,  since  Green's  theorem  can  be  used 
to  shift  one  derivative  from  u  to  v.  A  function  u  satisfying  (2.1.4)  is 
called  a  weak  solution  of  (2.1.3),  and  (2.1.4)  is  called  the  weak  form 
of  the  equation.  The  advantage  here  is  that  the  weak  solution  exists 
quite  generally  and  agrees  with  the  classical  solution  satisfying 
equation  (2.1.3)  pointwise  when  the  latter  exists. 
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The  condition  for  existence  and  uniqueness  of  the  solution  of 
(2.1.4)  is  that  the  problem  be  elliptic.  This  means  that  the  energy 
norm  |||  •  |||  defined  by 

ll|u|||  2  =  a(u,u) 

for  u  e  is  equivalent  to  the  first  Sobolev  norm  ||  •  ||  .  That  is. 


i  a2\\ 


u  e 


hJ(H) 

E 


for  fixed  cr^,  >  0  independent  of  u. 

If  the  problem  (2.1.3)  is  elliptic,  a  simple  calculation  shows 
that  for  £  a  non-negative  integer 


lull£+2  1  c  IILu!I£ 


U  e  H£+2(ft)  , 


for  some  c  >  0.  The  principal  results  in  elliptic  theory  go  the  other 
3  £ 

way.  For  £  >  -  -j  and  any  f  e  H  (fl)  the  weak  solutions  of  (2.1.1)  or 
(2.1.2)  satisfy  the  regularity  estimate 
(2-1.5)  l|u||£+2  5  c 


for  some  c  >  0,  assuming  that  the  domain  is  smooth  and  the  coefficients 

CO 

of  L  are  C  .  For  the  technical  definition  of  smooth  domains  see 

Babushka  and  Aziz  (1972)  or  Oden  and  Reddy  (1976). 

The  regularity  estimate  (2.1.5)  implies  existence,  uniqueness, 

and  that  the  solution  has  two  more  derivatives  in  the  mean  square  sense 

than  its  data.  In  particular,  if  f  is  infinitely  differentiable,  u  will 

be  as  well.  These  results  hold  only  for  smooth  domains.  For  domains 

that  are  only  piecewise  smooth  the  corners  introduce  singularities  that 

limit  the  regularity  of  the  problem.  On  a  convex  polygonal  domain  the 

2  £ 

solution  will  lie  in  H  but  not  necessarily  in  H  for  £  >  2.  In  this 


case  we  have  the  regularity  estimate 


t 


9 


14 

(2.1-6)  !Mla+2±  c(|f||A 

3 

for  -  —  <  H  <_  0.  For  non-convex  polygonal  domains  (2.1.6)  will  not  hold 

for  1  =  0  because  the  reentrant  comers,  that  is  corners  with  interior 

angles  exceeding  tt,  introduce  severe  singularities.  In  the  worst  case  of 

3  1 

a  slit  domain,  H  must  be  restricted  to  (-  -j,  -  -j) .  The  solution  has 
less  than  one  and  a  half  derivatives  in  the  mean  square  sense. 

2.2.  Finite  Element  Spaces 

The  finite  element  discretization  of  (2.1.3)  begins  with  the 
selection  of  a  finite  dimensional  subspace  M  of  Hj,(fi).  The  discrete 

Ej 

solution  to  (2.1.3)  is  taken  as  that  function  u  6  M  satisfying  the  finite 
element  equation 

(2.2.1)  a(u,<fi)  =  (f ,tp),  <p  e  M  . 

By  the  Ritz  theorem,  holding  for  self-adjoint  problems,  u  is  the  function 
in  M  nearest  to  the  true  solution  u  in  energy  norm. 

Finite  element  spaces  are  ordinarily  constructed  by  dividing 
the  domain  into  a  finite  family  of  elements,  usually  triangles  or 
rectangles.  Once  the  elements  are  chosen  M  may  be  taken  as  a  finite 
dimensional  space  of  piecewise  polynomial  functions  that  are  polynomial 
on  each  element.  Though  it  is  not  strictly  necessary  for  multi-level 
methods,  throughout  we  will  require  that  M  be  conforming,  that  is  that 
it  lie  in  H^,(fi).  This  translates  into  the  requirement  that  the  functions 
in  M  be  continuous  and  satisfy  the  essential  boundary  condition  (2.1.3b) 
exactly.  We  also  require  that  M  be  of  degree  at  least  one  in  order  to 
obtain  convergence.  For  a  finite  element  space  of  piecewise  polynomial 
functions  the  degree  is  the  integer  k  such  that  any  polynomial  of  degree 


ft 
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k  can  be  exactly  represented  on  every  element  away  from  the  boundary. 

Finally  it  is  necessary  to  bound  the  highest  degree  polynomial  on  any 

element  and  require  that  M  have  a  local  basis.  This  is  in  order  that  the 

basis  be  uniformly  linearly  independent.  The  question  of  the  linear 

Independence  of  the  basis,  or  conditioning  of  the  Gram  matrix  will  be 

taken  up  in  the  next  section. 

Approximation  questions  are  dealt  with  in  finite  element 

theory  in  terms  of  a  sequence  of  finite  element  spaces  of  increasing 

dimension  {M.}.  , ,  which  eventually  become  dense  in  the  sense  that  for 
1  1=1 

any  neighborhood  N  in  there  exists  J  such  that 

£j 

N  n  M,  ,1  o  for  j  >  J  . 

The  approximation  properties  of  each  of  these  spaces  M  depend  on  the 

size  of  the  elements  {e}  of  and  the  element  geometry,  as  well  as 

the  degree  of  the  finite  element  space.  Element  geometry  matters  only 

to  the  extent  that  it  must  not  degenerate  as  j  increases.  For  standard 

finite  elements  it  suffices  that  all  elements  in  M.  be  convex  and  have 

J 

interior  angles  bounded  away  from  0  and  it  uniformly  in  j. 


The  size  of  an  element  e  is  measured  by  its  diameter,  dg. 


defined  as 

d  =  sup  !|x-y||  , 

6  x,yee  JT 

where  ||  •  ||  „  denotes  the  Euclidean  norm.  The  family  of  spaces  (M  l  is 
l  2 

said  to  be  quasi-uniform  if  there  is  a  corresponding  family  of  parameters 

{ hj }  and  a  >  1  such  that 


(2,2.2) 


■r  h  <  d  <  o  h. 
a  j  -  e  -  j 


for  all  elements  e  of  Mj  and  all  j.  If  not,  that  is  if  the  constant  0 
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must  depend  on  j,  the  family  of  spaces  {M ^ }  is  said  to  be  locally  refined. 
While  strictly  speaking  these  concepts  apply  only  to  an  infinite  family 
of  finite  element  spaces,  we  will  talk  about  an  individual  finite  element 
space  being  quasi-uniform  or  locally  refined.  This  usage  is  extremely 
convenient  though  it  should  be  born  in  mind  that  this  really  makes  sense 
only  with  respect  to  membership  in  an  infinite  family  of  spaces {Mj}. 

This  definition  of  the  mesh  parameter  h^  associated  with  a 
quasi-uniform  finite  element  space  M  is  not  completely  standard.  Often 
hj  is  taken  as  the  maximum  of  the  element  diameters.  However,  the 
extra  freedom  allowed  by  (2.2.2)  makes  no  difference  in  the  error 
estimates  since  a  is  a  fixed  constant,  and  is  quite  convenient.  In 
particular  this  freedom  makes  it  easy  to  have  a  quasi-uniform  family 

1)  Parameter:Lzed  by  h  rather  than  j.  In  this  case  the  diameters 
of  elements  in  are  required  to  satisfy 

^  h  <  d  <  oh  . 

o  —  e  — 

The  standard  finite  element  error  estimates  give  the  rate 
of  convergence  in  any  Sobolev  norm  of  the  finite  element  approximations 
in  a  quasi-uniform  space  M  in  terms  of  the  mesh  parameter  h  and  the 
degree  of  M^.  For  smooth  solutions  the  convergence  rate  may  be  made 
arbitrarily  high  by  taking  finite  element  spaces  {M^}  of  sufficiently 
high  degree.  But  on  a  polygonal  domain  the  corner  singularities  limit 
the  rate  of  convergence  attainable  by  any  family  of  piecewise 
polynomial  finite  element  spaces. 

The  most  important  error  estimate  in  the  next  chapter  will  be 


the  L  error  estimate 


« 
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(2.2.3) 


I  u-uh||  <  Ch2  ||f|| 
L  L 


which  holds  for  quasi-uniform  spaces  {Af  }  of  degree  at  least  one  when 

h 

p 

the  partial  differential  equation  is  H'  regular,  that  is,  when  (2.1.6) 
holds  with  S.  =  0.  Thus,  this  bound  applies  to  the  Dirichlet  problem 
(2.1.1)  and  the  Neumann  problem  (2.1.2)  on  any  convex  polygonal  domain. 

The  convergence  in  energy  for  these  two  problems  is  of  first 
order  on  any  convex  polygonal  domain.  On  non-convex  domains  it  is  of 
less  than  first  order.  For  a  domain  with  a  slit,  the  worst  case 
ordinarily  arising,  it  is  of  order  one-half. 

Much  of  the  utility  of  the  finite  element  method  comes  from 
the  fact  that  these  meager  convergence  rates  attained  on  polygonal 
domains  are  easily  improved  through  the  addition  of  singular  functions 
or  by  doing  local  refinement  to  cope  with  the  corner  singularities.  In 
chapter  four  the  use  of  multi-level  methods  on  locally  refined  finite 
element  grids  will  be  considered.  A  simple  error  estimate  showing 
the  quality  of  approximation  of  finite  element  methods  on  such  grids 
will  be  given  there. 


2.3.  Linear  Equations 

By  selecting  a  basis  for  the  finite  element  space  Af,  the 
finite  element  equation  (2.2.1)  is  translated  into  a  linear  system. 

Let  Af  be  a  finite  element  space  for  problem  (2.1.3)  and  let 

be  a  basis  for  Af.  Any  function  u  e  Af  can  be  expressed  in  this  basis  as 

N 

u  =  Z  u  y.  , 
i=l 
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where  ju  is  the  coordinate  vector  of  u  relative  to  this  basis.  The 
convention  of  using  underbars  for  vectors  will  be  adhered  to  throughout. 

Given  the  basis  we  can  define  the  mass  and  stiffness 

matrices  M  and  K  with  elements 

"«  '  <V  V 
kU  ‘  ,<*i’  V 

respectively.  M  is  just  the  Gram  matrix  of  the  basis  relative  to  the 

2 

L  inner  product  while  K  is  the  Gram  matrix  in  the  energy  inner 
product.  To  form  the  linear  system  corresponding  to  the  Ritz  equation 

(2.2.1) ,  let  f  e  L2(Q)  and  let  £  be  given  by 

If  =  (f,  4^),  i  -  1 . N  • 

Then  (2.2.1)  is  equivalent  to  the  linear  system 

(2.3.1)  Ku  =  F 

where  ii  is  the  coordinate  vector  of  the  finite  element  solution. 

Numerical  solution  of  the  linear  system  (2.3.1)  is  aided  by 
several  properties  the  matrices  K  and  M  have.  First,  M  will  always  be 
symmetric  positive  definite  and  since  L  is  self-adjoint  K  will  be 
symmetric  positive  definite  as  well.  Second,  for  the  usual  finite 
element  bases  M  and  K  will  be  quite  sparse.  Finally,  both  M  and  K  are 
usually  well  conditioned.  The  common  finite  element  bases  are  largely 
motivated  by  these  last  two  considerations  which  are  of  great  importance 
for  the  numerical  inversion  of  (2.3.1)  or  (2.3.2)  whether  this  is  done 
directly  or  iteratively. 

To  make  these  considerations  precise,  consider  a  family  of 

finite  element  spaces  {M,}  of  increasing  dimension.  We  require  that 

J  J"1  V 

(i )  i 

each  space  M  ^ ,  1  _<  j  <  00 ,  have  a  basis  that  is  local  in  the 


« 
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sense  that  the  number  of  non-zero  elements  in  each  row  of  the  corresponding 
mass  and  stiffness  matrices  is  bounded  uniformly  in  j.  This  will  be  true 

C)  N' 

if  the  number  of  basis  functions  in  {<t>.  t."*,  whose  supports  intersect 

1  i=l 

the  supports  of  any  given  basis  element  <j>^^  is  bounded  uniformly  in  j. 

Usually,  under  these  conditions  the  condition  number  k(M.)  of 

3 

the  mass  matrix  ,  1  _<  j  <  °°,  will  be  uniformly  bounded.  For  this  to 
hold  the  proper  scaling  must  be  chosen  and  the  condition  numbers  of 
the  element  mass  matrices  must  be  uniformly  bounded.  This  is  true  for 
the  common  finite  element  bases  as  long  as  the  element  geometry  does 
not  degenerate.  Then  if  the  basis  elements  are  scaled  as 

ll<t>p)H  =1.  1  1  i  £  N  ,  1  1  j  <  ^  , 

which  we  assume  throughout,  it  is  immediate  that  the  eigenvalues  of 
IF  will  lie  in  an  interval  [^,  a]  for  a  >  1  independent  of  j. 

Since  the  stiffness  matrices  {K.}  are  exactly  analogous  to 
the  matrices  arising  in  finite  difference  discretizations,  their 
condition  numbers  cannot  be  uniformly  bounded.  Just  as  in  the  finite 
difference  case  the  maximum  eigenvalue  of  K.  satisfies 


X  (K.)  <  ch  . 
max  j  —  min 


(2.3.2) 

under  ordinary  circumstances,  where  h 


min 


is  the  diameter  of  the 


smallest  element  of  .  The  minimum  eigenvalue  of  can  be  bounded 
below  by  a  simple  argument.  We  have 


T  T  T 

x  K.x  x  K  x  x  M  x 

A^K  )  =  min  — >  min  ■  --  ^  min  — 

XX  x  Mj2c  XX 


A1(m"  K.)  AX(M  ) 


>  ^(Mj1  K  )/0  . 
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The  operator  M  ^  is  the  linear  operator  corresponding  to  the  finite 
element  equation  (2.2.1).  Since  the  finite  element  eigenvalues 
approximate  the  eigenvalues  of  the  PDE  from  above,  X^(K^)  must  be 
bounded  below  uniformly  in  j.  Thus  one  ordinarily  has 


(2.3.3) 


k(K.)  <  ch 
J  ~ 


-2 

min 


or  for  a  quasi-uniform  spaces 

(2.3.4)  kO^)  <  ch' 


The  difficulty  in  iteratively  inverting  the  linear  system  (2.3.1)  arises 
directly  from  this  growth  of  the  condition  of  the  stiffness  matrix  at, 
h  decreases. 

These  bounds  on  the  condition  of  M  and  K  are  proved  in 
Fried  (1971)  and  may  be  found  also  in  Strang  and  Fix  (1973). 


21 


3.  QUASI-UNIFORM  GRIDS 

In  this  chapter  a  multi-level  algorithm  simpler  than  those 
considered  in  previous  theoretical  work  will  be  given.  It  will  be 
shown  to  be  of  optimal  order,  producing  a  good  solution  on  a 
quasi-uniform  finite  element  grid  in  0(N)  operations  on  an  N 
point  grid. 

3.1.  Introduction 

In  this  chapter  we  look  at  a  multi-level  algorithm  for 
self-adjoint  elliptic  boundary  value  problems  on  quasi-uniform  finite 
element  grids.  This  algorithm  is  related  to  algorithms  considered  by 
Nicol aides.  Bank  and  Dupont,  Hackbusch  and  others,  but  is  less  complex 
than  these  other  algorithms.  Like  many  of  those  considered  by  others, 
the  algorithm  here  is  of  optimal  order,  producing  a  second  order 
accurate  solution  in  0(N)  operations  on  an  N  point  grid. 

The  optimal  order  convergence  of  the  algorithm  in  this  chapter 
was  originally  shown  for  finite  difference  discretizations  by  a 
cumbersome  Fourier  technique,  and  then  for  finite  element  discretizations 
in  the  same  way.  Here  an  easier  proof,  following  the  ideas  of  the  above 
authors,  will  be  given.  The  proof  is  based  particularly  on  the  work  of 
Bank  and  Dupont  but  employs  several  minor  simplifications.  Throughout 
we  follow,  for  the  most  part,  their  notation.  It  is  convenient  here 
since  little  reference  to  algebraic  quantities  such  as  vectors  or 


matrices  will  be  made. 


In  this  section  and  the  next  we  examine  the  convergence  of  a 
simple  multi-level  algorithm.  This  algorithm  can  be  used  to  obtain 
finite  element  approximations  to  the  solutions  of  the  elliptic 
problems  (2.1.1)  and  (2.1.2)  on  a  convex  polygonal  domain  Q.  The 

restriction  to  convex  domains  is  necessary  since  we  need  to  use  the 

2  2 
L  error  estimate  (2.2.3)  which  requires  H  regularity.  Using  curved 

elements  as  in  Bank  and  Dupont  [1978]  the  convergence  results  given 

here  extend  immediately  to  piecewise  smooth  domains  free  of  reentrant 

corners.  For  details  on  the  treatment  of  curved  elements,  the 

interested  reader  is  referred  to  their  paper. 


The  algorithm  considered  here  uses  a  nested  family  of 

CO 

quasi-uniform  finite  element  spaces  ,  nested  in  the  sense  that 

J  J 

M,  c  M  ,  1  <  i  <  00  . 

j  1+1 


To  approximately  solve  the  finite  element  equations  for  one  of  these 
spaces  k  _>  2,  this  algorithm  uses  a  simple  simultaneous  displacement 
smoothing  iteration  applied  to  the  equations  for  plus  approximate 
solution  of  related  problems  in  If  k  -  1  >  2  these  approximate 

solutions  are  obtained  recursively  by  applying  the  same  algorithm  to 
the  equations  for  The  recursion  continues,  finally  requiring 

approximate  solutions  to  the  equations  for  These  can  be  found 
directly  or  by  some  convenient  iteration. 


The  multi-level  method  considered  here  is  related  to  algorithms 


for  finite  elements  by  Nicolaides,  Hackbusch,  and  Bank  and  Dupont  and  to 
methods  for  finite  difference  equations  considered  by  Brandt,  Hackbusch, 
Nicolaides,  Bakhvalov,  Fedorenko  and  others.  The  algorithm  here  differs 


from  algorithms  described  by  others  in  that  the  solution  to  the  equations 

for  requires  only  two  approximate  solutions  to  the  equations  for 

In  other  multi-level  algorithms,  to  solve  the  equations  for 

approximate  solution  of  related  problems  in  is  used  in  two  ways. 

First  it  is  used  as  a  means  of  obtaining  a  good  starting  value  for 

iteration  on  Subsequently  it  is  used  as  an  acceleration  procedure 

for  iteration  on  M,  . 

k 

The  first  approximate  solution  of  the  equations  for 
is  used  to  obtain  a  starting  value  for  iteration  on  and  considerable 
accuracy  is  required.  It  is  necessary  to  reduce  the  convergence  error 
of  the  problem  on  ^  so  it  is  comparable  to  the  truncation  error  on 
j .  in  subsequent  solutions  of  related  problems  in  ^  used  to 
accelerate  convergence  of  the  iteration  in  less  accuracy  is  required. 
For  example,  in  one  of  the  proofs  given  by  Bank  and  Dupont  and  in  the 
related  proof  of  theorem  4.1  in  chapter  four  of  this  thesis  it  suffices 
to  solve  the  related  problem  with  a  relative  accuracy  of  33%. 

The  point  of  this  chapter  is  that  if  all  of  the  related 
problems  in  M  ^  are  solved  to  the  same  accuracy,  namely  the  error  in 
the  approximate  solutions  must  be  comparable  to  the  truncation  error 
in  M  then  only  two  approximate  solutions  in  are  required. 

To  approximately  solve  the  elliptic  problem  (2.1.1)  or  (2.1.2)  in 
M^,  k  >_  2,  the  algorithm  begins  by  obtaining  a  good  approximate 
solution  to  the  finite  element  equations  for  This  approximate 

solution  is  used  as  a  starting  value  for  a  smoothing  iteration  on 
Though  this  iteration  could  be  continued  until  adequate  convergence 
was  attained  this  would  be  inefficient.  Instead,  after  a  fixed  number 
of  iterations,  m,  independent  of  k,  we  terminate  the  iteration. 


The  resulting  approximate  solution  is  then  corrected  by  subtracting  from 
it  the  projection  of  the  convergence  error  onto  This  corrected 

approximate  solution  in  is  then  taken  as  the  accepted  solution. 

The  reason  this  is  effective  goes  back  to  the  work  of 
Southwell.  He  observed  that  at  least  initially  relaxation  rapidly 
reduces  the  residual,  while  the  error  itself  may  be  only  slowly  reduced. 

At  the  end  of  our  algorithm's  smoothing  iteration  the  residual  will  be 
small  relative  to  the  convergence  error.  This  means  that  the 
convergence  error  is  quite  smooth.  In  this  case  the  convergence  error 
can  be  well  approximated  in  the  coarser  space  ^  and  the  final 
correction  will  remove  most  of  it. 

Brandt  describes  this  same  phenomenon  in  terms  of  the  reduction 
of  Fourier  components  of  the  convergence  error.  Most  convergence 
proofs,  including  ours,  follow  Fedorenko  who  considered  the  reduction 
of  the  error  components  in  the  direction  of  eigenfunctions  of  the 
discrete  partial  differential  operator.  Fedorenko  refers  to  "good" 
meaning  smooth  eigenfunctions  and  "bad"  meaning  oscillatory  eigenfunctions. 

For  practical  purposes  many  possible  smoothing  iterations 
would  be  acceptable.  Jacobi,  Gauss-Seidel ,  and  symmetric  SOR  are  most 
often  used.  The  author  expects  that  improvement  could  be  made  through 
the  use  of  Chebyshev  acceleration,  although  in  numerical  tests  so  far 
changing  to  more  subtle  iterative  schemes  has  had  little  effect  for 
the  well  behaved  self-adjoint  problems  considered  here  (Bank  and 
Sherman  [1979],  Brandt  [1977]). 

In  this  section  the  smoothing  iteration  will  be  taken  as  a 
simple  but  computationally  unattractive  simultaneous  displacement 
iteration  considered  also  by  Bank  and  Dupont.  The  advantage  of  this 


t 


iteration  is  that  the  resulting  convergence  is  easy  to  analyze  and  is 
basis  independent.  In  the  next  section  we  will  show  that  this  smoothing 
iteration  can  be  replaced  by  a  more  computationally  attractive 
iteration  without  altering  the  convergence  result.  In  this  case  the 
algorithm  will  be  shown  to  be  of  optimal  order,  producing  a  second 
order  accurate  solution  in  0(N)  operations  on  a  finite  element  grid 
with  N  elements. 

bet  {M.}“  ,  be  a  nested  family  of  quasi-uniform  finite 
1  J  =  1 

element  spaces,  and  suppose  we  have  associated  mesh  parameters  {h. } 
satisfying 


0.2.1) 


h.  =  P1“jh1  ,  j  >  1  , 


for  p  >  1  a  fixed  constant.  Let  u^  €  M^. ,  j  _>  1,  be  corresponding 
finite  element  solutions  of  problem  (2.1.1)  or  (2.1.2).  That  is, 
u('>)  E  Mj  are  functions  satisfying 


(3.2.2) 


a(uVJ\<f>)  =  (f,$)  ,  <f>  €  M 


From  the  error  estimate  (2.2.3) 


(3.2.3) 


|u-„(3)||£c  ||f  I 


for  fixed  c  >  0  where  u  is  the  true  solution.  In  the  space  associated 

with  the  finite  element  equation  (3.2.2)  we  have  eigenfunctions 
C)  N. 

(ill.  0.~'  ,  where  N.  is  the  dimension  of  M.,  and  associated  eigenvalues 
i  i=l  ]  J 

...  N. 

{Y  }i=i wlth 

0  <  A«>  <  A<j)  <...<  \<j)  . 

1  —  L  —  —  N . 

J 


That  is 


« 
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aO^.flO  =  ,  <J)  e  , 

for  1  £  i  £  N  .  From  the  eigenvalue  bound  (2.3.2)  and  the  fact  that 
the  eigenvalues  of  1L  ,  the  mass  matrix  corresponding  to  ,  lie  in 
[^,  a]  for  a  >  1  independent  of  j,  we  have  the  bound 

(3.2.4)  A„  £  c  hT2  . 

j  2 

To  solve  the  finite  element  equation  in  k  £  2,  the 

algorithm  to  be  considered  produces  a  sequence  of  iterates  u^, 

u^^  in  beginning  with  an  approximate  solution  u^  to  the  finite 

element  equations  in  For  some  fixed  positive  integer  m,  the 

approximate  solution  u  will  be  taken  as  the  final  answer.  The 

ntTi 

iterates  u^,  1  <  £  _<  m,  are  produced  by  a  simple  smoothing  iteration 

while  the  last  one,  which  will  be  the  accepted  solution,  is  the  result 

of  a  correction  involving  the  approximate  solution  of  a  related 

2 

problem  in  For  a  fixed  constant  c^  >  0  and  data  f  in  L  (fi)  the 

algorithm  is  as  follows: 

1.  Let  Uq  G  ^  be  an  approximate  solution  to  (3.2.2) 

satisfying 


(3.2.5) 


lvu(k_1)ll  -  co  (hk-i)2  11  f| 


Ck-ll 

where  uv  is  the  solution  of  (3.2.2), 


2.  Let  u? ,  1  <  £  £  m  be  the  element  of  M  satisfying 


(3.2.6)  (ujl-u^_1,4i) 


(a(ug,_i*<fr)  “  * 


3.  Let  v  £  M,  ,  be  the  solution  of 
k-1 


4)  e  Mk  . 


a(v,<J>)  -  a(um,4>)  -  (f,<J>)  ,  <f>  6  Mfc_1  . 


9 


That  is ,  v  is  given  by 


i(v,<J>)  =  - (r ,<{>)  ,  <j>  e  M 


where  r  C  f.l  is  defined  by 
k 


=  (f ,<f>)  -  a(u  ,4>)  , 


e  M  . 
k 


Let  v  be  an  approximation  to  v  satisfying 


(3.2.7) 


|v-v||  <  cQ  hk_x  II  r  I 


and  let  u  be  given  by 
m+  L 


(3.2.8) 


u  ,  ,  =  u  -  V 
m+1  m 


The  corrected  approximate  solution  u  ^  is  taken  as  the  final  answer. 

If  we  can  show  that  the  approximate  solution,  u^^,  produced 

is  accurate  enough  uniformly  in  k  _>  2,  the  use  of  this  algorithm 

recursively  to  solve  the  subproblems  in  steps  one  and  three  will  be 

immediately  justified.  This  is  the  content  of  the  following  theorem. 

Theorem  3.1.  For  any  c^  >  0  there  is  a  positive  integer  m, 

independent  of  k,  such  that  the  approximate  solution  u  . ,  produced  by 

m+1 

algorithm  ( 3. 2 . 5)- (3. 2 . 8)  satisfies 


(3.2.9) 


I Untfl  '  u(k>  H  <  CQ  hk  !ifll  ’ 


Proof.  For  convenience  we  drop  the  superscript  k  on 


eigenfunctions  and  eigenvalues  and  the  dimension  of  M^.  From  (3.2.3) 
and  the  triangle  inequality 

![„<»  -  u<k'1)||<  c<h2  *  hjhp  ||f|| 

<  c(l  +  c2>  h2  ||f  II  . 


f 


« 
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From  this  and  (3.2.5) 

(3.2.10)  ||u0  -  u(k)  ||  <  (c(l  +  p2)  +  cQ  p2)  h2  ||f||  . 

Now,  for  SL  =  0,  1,...,  ra+1  let  e^  be  the  convergence  error 


Expanding  e^  in  eigenfunctions 


e 

0 


N 

Z 

i=l 


CL.  <p. 

X  1 


Then  by  (3.2.6) 


e 

m 


N  A. 

*  V1  -  4>m 

x=l  N 


It  remains  to  estimate  the  result  of  the  correction  In  step  3. 


We  have 


N  A. 

r  =  Z  a.  A  (1  -  y±)m  ip 
i=l  11  AN  i 


and  by  the  error  estimate  (2.2.3) 


||v  -  em||  <  c  h2_1  || r ||  ; 

so  by  (3.2.7) 

1 1  1 1  =  llem-vll  £  llem-v||  +  !lv-vll 

£ch2_i  ||  r  ||  +  cQ  h2^  ||  r  || 

<  (c  +  cQ)  P2  h2  II  r II 

«  (c  +  cQ)  p2  h2  ||^  ot  A±(l  -  y^)m  ^|| 

<  (c  +  c.)  p2  h2  II e_ II  max  (A (1  -  y~)m) 

0  0<A<Am  an 

- N 
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where  the  last  step  depends  on  the  L  orthogonality  of  the  eigenfunctions 

{li'j).  By  a  simple  calculation 

/  \  f  X  ,m.  ^  N 

max  (A  (1  -  — )  )  <  — —  . 

0<A<,\  AN  nH'1 

—  —  N 


Then  from  the  eigenvalue  bound  (3.2.4) 


iemhlll-  C  p2  "1ST2  He0l 


It  follows  from  this  and  (3.2.10)  that 


I  Villi  c  o2  <C<1  +  D2)  +  c0  p2)  ||f  II  . 


Choose  m  large  enough  that 

2 

(3.2.11)  m+1  >  ^-2-  (c  +  c  )  (c(l  +  p2)  +  c  p2)  , 

0  U 

and  the  result  follows.  □ 

In  this  proof  we  have  carried  the  constants  p  and  c^  explicitly 

to  show  the  dependence  of  m  on  these  quantities.  From  (3.2.11)  m 
4 

apparently  goes  as  —  .  This  sharp  dependence  of  the  number  of  iterations 
C0 

m  on  the  mesh  ratio  is  typical  of  multi-level  algorithms,  but  the 
dependence  on  c^  is  not.  For  more  common  multi-level  algorithms,  such 
as  that  considered  by  Nicolaides  [1977]  and  that  considered  in  chapter  4 
for  locally  refined  grids,  m  depends  on  c^  only  as  |log(Cg)|.  Thus  if 
one  wishes  to  reduce  the  convergence  error  well  below  the  truncation 
error,  the  algorithm  here  cannot  be  recommended.  However,  in  the  usual 
case  where  it  suffices  to  have  convergence  error  comparable  to  the 
truncation  error  this  algorithm  should  perform  about  as  well  as  other 


multi-level  methods. 
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3.3.  Computational  Cost 

While  iteration  (3.2.6)  used  in  step  two  of  algorithm  (3.2.5)- 
(3.2.8)  was  convenient  in  the  last  section  since  it  made  the  analysis 
there  basis  independent,  it  is  impractical.  Each  of  the  iterations, 
(3.2.6),  required  the  solution  of  a  linear  system  involving  the  mass 
matrix  because  of  the  implicit  definition  on  the  left  side.  To  avoid 
this,  we  will  show  in  this  section  that  the  smoothing  iteration  (3.2.6) 
can  be  replaced  by  a  simple  simultaneous  displacement  iteration  without 
impairing  the  convergence  result.  This  analysis  follows  Bank  and 
Dupont  [1978].  With  this  replacement  the  overall  algorithm  is  of  optimal 
order,  producing  a  second  order  accurate  solution  in  0(N)  operations  for 
a  finite  element  grid  with  N  elements. 

Given  a  basis  for  M^,  for  each  u  6  there  corresponds 

N 

a  coordinate  vector  u  6  IR  .  As  in  the  last  section,  we  have  temporarily 
dropped  the  subscript  k  on  the  dimension  N  of  and  on  eigenvectors  and 
eigenvalues.  Let  b(*,  •)  denote  the  bilinear  form 


T 

b(u,v)  =  u  v 

The  norm  induced  by  this  bilinear  form  is  the  same  as  the  Euclidean  norm 
on  corresponding  vectors, 

l|u|£  =  b(u,u)  =  !|u||22  . 

Now  let  {^}  be  the  eigenvectors  of  the  stiffness  matrix  K  with 
corresponding  eigenvalues 

0  <  X  <  <. . .<  I  . 

i  —  l  —  —  N 

The  eigenvectors  {j^}  are  coordinate  vectors  for  eigenfunctions  {^} 


satisfying 
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(3.3.1) 


a(iK,  4>)  =  Ai  bC^,  <p)  , 


and  the  orthogonality  relations 

a(ijk,  ifO  =  b^,  $  )  =  0  for  i  ?  j  . 

We  wish  to  consider  algorithm  (3. 2 . 5)- (3. 2. 8)  of  the  last 
section  but  with  the  smoothing  iteration  (3.2.6)  replaced  by  the 
simultaneous  displacement  iteration 


(3.3.2)  b(u„  -  u0  1  ,  <p)  =  -  —  (a(uJl_1,  <f>)  -  (f,  <P))  , 


l-l 


'N 


e  M, 


The  use  of  A^  here  is  only  for  convenience.  It  could  be 
replaced  by  a  convenient  upperbound  on  the  eigenvalues  of  the  stiffness 
matrix,  such  as  (2.3.2),  without  altering  the  convergence  result. 

The  convergence  of  this  algorithm  is  shown  by 

Theorem  3.2.  For  any  c^  >  0  there  is  a  positive  integer 


m  independent  of  k  >  2  such  that  the  approximate  solution  u 


(k) 

nrt-1 


produced  by  algorithms  (3. 2. 5 (- (3. 2 . 8)  with  (3.2.6)  replaced  by  (3.2.2) 
sat  is f i es 


m+1 


“(k)Hi  hj  II  f  I 


Proof.  The  proof  is  a  slight  modification  of  the  proof  of 
Theorem  3.1.  Expanding  in  the  eigenfunctions  instead  of  in  the 


eigenfunctions  {ii^}  we  have 


Then,  as  before, 


e  =  E  a  ♦.  . 
U  i=l 


N  A . 

E  a  (1  -  - — )m  . 

i=l  1  ANk 


'4 


f 
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Let  r€  M  be  given  by 


(r,<p)  =  -  a (e  ,  <j>)  ,  <t>  e  M 

m  it 


and  r  G  M  by 


Then 


so  we  have 


b(r ,$)  =  -  a(e  ,  <j>)  , 
m 


r  =  Mr  ; 


e  Mk  . 


2  T  —  1  T  _  "j 

I  r  ||  =  r  Mr  =  (M  ?)  M(M  ?) 


~T  M-l. 
=  r  M  r 


<o||tl£2-o||r|g 


As  in  the  proof  of  Theorem  3.1 


i enH.1  II  i  (cd  +  p2)  +  c0  p2)  IMI 

<  (C (1  +  P2)  +  C0  p 2)01/2  hk  llrllh 


Now 


N 


i.m 


'H.  -II  L  “i  V1  -  tr  *i  "b 

i=l  AN 


<  || e  |L  max  (A(l  - 

0  ^  0<X<\ 

- -  N 


where  the  first  inequality  uses  orthogonality  of  the  eigenfunctions  {^} 
in  the  inner  product  b(*,  •).  Then  since 


* 


(3.3.3) 


T  M  _  ,m1/2  T  .1/2  . 

«0  »«0  -  <«  <M  %) 


-ll»1/2i0IP2ii||e0!P2-i||.0l!2b 

JO 


and  since  the  maximum  eigenvalue  satisfies  the  bound  (2.3.2)  the  rest 
of  the  proof  goes  through  exactly  as  before.  □ 

This  theorem  shows  that  by  twice  approximately  solving  the  lower 
dimensional  linear  system  for  and  doing  a  fixed  number  of  simultaneous 

displacement  iterations  on  one  can  obtain  a  second  order  accurate 
solution  in  M  .  Since  the  convergence  shown  by  Theorem  3.2  is  independent 
of  k  j>  2,  the  two  related  problems  in  ^  can  be  solved  by  recursively 
applying  the  same  algorithm  to  their  smaller  linear  systems.  Continuing 
the  recursion  one  ends  up  repeatedly  having  to  solve  the  linear  system  for 
M^.  This  can  be  done  either  directly  or  by  any  convenient  iterative  method. 

To  estimate  the  computational  cost  of  this  algorithm  observe 


that  from  (3.2.1)  the  dimensions  of  the  finite  element  spaces  (Mj }  must 
satisfy 


:l  P2^  i  Nj  1  c2  p2''  > 


11  j  < 


for  c^,  >  0  independent  of  j.  This  assumes  that  h^,  the  mesh  size  of 

is  fixed  and  the  degree  of  the  highest  degree  polynomial  on  any  element 
is  bounded  uniformly  in  j.  Then,  if  s^  is  the  storage  required  for  the 
direct  or  iterative  inversion  of  the  linear  system  for  the  total 
storage  must  be  proportional  to 


S0  +  *  Nj  -  S0  +  C2  *  P 

J-2  J  j=2 


5  s° +  TP 


<— *r2>»k 

u  1  1-p  k 


for  p  >  1.  Since  s^  is  fixed,  assuming  the  mesh  ratio  p  is  greater  than 

one,  the  total  storage  required  depends  only  linearly  on  the  dimension  of 

the  solution  space  M,  . 

k 

To  bound  the  computation  time,  observe  that  the  time  required 
to  do  one  of  the  simultaneous  displacement  iterations  (3.3.2)  for  M.  is 
proportional  to  N  .  Let  t^  be  the  computation  time  for  the  approximate 
solution  of  the  linear  system  for  M^.  Then  taking  into  account  that  we 
must  recursively  apply  the  algorithm  twice  to  MR  ^ ,  four  times  to  M,  ^ , 
and  so  on  the  computation  time  T  is  bounded  by 


(3. 3.4) 


1  *v  1  • 

T  <  2  t.  +  cm  l  2  3  N. 

-  °  j=2  3 


<  2k  t  +  cc  m  Z  2k~j  p2k 

3=2 

<  2k  tQ  +  cc2  mp2k  Z  (-^)k"j  . 

3=2  P 

We  have  neglected  the  computational  cost  of  carrying  out  the 
mappings  of  functions  on  one  grid  to  another  in  steps  one  and  three  of 
algorithm  ('3.2.  5)— (3.2.8).  Despite  the  fact  that  the  spaces  }  are 
nested,  the  mappings  required  in  steps  one  and  three  are  not  free  since 
in  general  each  space  M  uses  a  different  basis.  It  is  however  easy  to 
see  that  the  above  bound  on  the  computation  time  is  unaltered  by  taking 
into  account  this  additional  cost. 

Now  from  (3.3.4)  we  have  the  bounds 


0(Nk> 


for  p  >  /2 


0(Nr  log (NR) )  for  p  =  /2 


log  2 

2  log  p 


for  1  <  p  <  /2  . 
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The  natural  choice  p  =  2  gives  an  optimal  order  algorithm  and  is  convenient 
in  programming. 


4.  LOCALLY  REFINED  GRIDS 

This  chapter,  which  contains  the  principal  results  of  this 
thesis,  presents  two  multi-level  algorithms  achieving  0(N)  complexity 
bounds  on  locally  refined  grids.  One  of  these  attains  this  optimal 
complexity  bound  under  weaker  than  expected  restrictions  on  the 
dimensions  of  the  finite  element  spaces  used  by  the  method.  These 
optimal  order  convergence  results  are  a  consequence  of  an  approximation 
theorem  proved  in  this  chapter.  Aside  from  its  use  in  extending  multi¬ 
level  convergence  theory  to  locally  refined  grids,  this  approximation 
result  is  of  interest  since  it  relies  only  on  local  properties  of  the 
finite  element  spaces  used  and  is  not  based  on  elliptic  regularity. 

As  a  consequence  it  provides  an  independent  verification  of  the  0(N) 
convergence  of  multi-level  methods  on  non-convex  domains  shown  previously 
by  Bank  and  Dupont. 

4 . ] .  Introduction 

There  are  a  variety  of  reasons  why  one  might  wish  to  use  locally 
refined  grids.  They  can  be  used  to  maintain  high  order  accuracy  in  the 
face  of  singularities  caused  by  discontinuous  coefficients  or  the  corners 
of  the  domain,  or  to  resolve  boundary  layers.  They  may  simply  result 
from  an  adaptive  discretization.  They  jean  also  be  used  to  give  good 
resolution  of  the  part  of  the  solution  of  greatest  interest  without 
incurring  the  computational  cost  of  a  fine  global  grid.  This  situation 
occurs,  for  examp’e,  on  unbounded  domains.  As  long  as  the  ratio  of  the 
diameter  of  the  largest  element  to  that  of  the  smallest  remains  uniformly 
bounded  on  the  family  of  finite  element  spaces  of  interest,  multi-level 
convergence  results  for  quasi-uniform  grids,  such  as  those  in  the  last 
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chapter,  remain  valid.  However,  in  many  cases  of  interest  this  ratio  of 
element  diameters  becomes  unbounded  on  the  family  of  finite  element 
spaces  being  used.  In  this  case,  while  results  for  quasi-uniform 
spaces,  such  as  those  in  the  last  chapter,  yield  an  0(N)  work  bound  for 
each  finite  element  space,  the  constant  hidden  in  the  0(N)  work  bound 
will  depend  on  the  ratio  of  element  diameters.  Such  bounds  are  clearly 
unsatisfactory  for  locally  refined  finite  element  spaces. 

While  nearly  all  of  the  convergence  results  in  the  literature 
apply  only  to  the  case  of  quasi-uniform  finite  difference  or  finite 
element  grids,  multi-level  methods  are  being  used  increasingly  for 
locally  refined  grids,  Brandt  (1977),  Bank  and  Sherman  (1978).  In  fact, 
one  of  their  principal  advantages  is  that  their  rate  of  convergence 
remains  very  good  on  locally  refined  grids.  By  contrast  other  iterative 
methods,  whose  rates  of  convergence  depend  on  the  condition  number  of 
the  stiffness  matrix,  do  poorly  on  locally  refined  grids.  For  example, 
the  number  of  iterations  required  with  conjugate  gradient  iteration 
goes  as  /K(K),  where  K( K)  is  the  spectral  condition  number  of  the  stiffness 
matrix,  Axelsson  (1977).  This  condition  number  satisfies 

K(K)  >  c  h”2 

—  min 

where  h  .  is  proportional  to  the  diameter  of  the  smallest  element  in  the 
min 

grid.  Thus,  X(K)  may  grow  much  faster  as  a  function  of  the  number  of 
unknowns  on  a  locally  refined  grid.  This  presents  a  very  severe 
difficulty  for  nearly  all  iterative  methods.  As  will  be  shown  in  this 
chapter,  multi-level  methods  do  not  suffer  from  this  difficulty. 

The  algorithms  in  this  chapter  differ  considerably  from  that  in 
the  last  chapter,  aside  from  their  application  to  locally  refined  grids. 

In  the  algorithm  of  the  last  chapter,  to  approximately  solve  the  elliptic 


I 


-I 


* 


» 


38 

problem  on  the  finest  grid,  a  related  problem  on  the  next  coarser  grid 
was  solved  twice.  It  was  solved  once  at  the  beginning  to  obtain  a 
starting  value  for  the  simultaneous  displacement  iteration  and  once  at 
the  end  as  a  final  correction.  In  the  algorithms  ofthis  chapter  and 
In  other  algorithms  in  the  literature,  solution  of  a  related  problem  on 
the  next  coarser  grid  is  used  more  frequently,  typically  after  every 
few  smoothing  iterations  on  the  finest  grid.  Thus,  the  basic  operation 
in  such  an  algorithm  is  best  thought  of  as  a  multi-level  iteration 
consisting  of  a  simple  relaxation  iteration  periodically  accelerated 
by  the  solution  of  a  related  problem  on  the  next  coarser  grid. 

Algorithms  like  those  considered  in  this  chapter  have  been 
described  for  quasi-uniform  grids  by  Bank  and  Dupont,  Nicolaides, 
Hackbusch,  and  others.  They  are  also  related  to  methods  for  locally 
refined  grids  suggested  by  Brandt  and  very  closely  to  the  method  used  in 
the  adaptive  finite  element  code  of  Bank  and  Sherman.  An  algorithm 
similar  to  that  of  chapter  three  could  have  been  used,  but  probably 
would  not  have  been  as  practical  as  the  algorithms  here.  One  could  argue 
that  in  any  case  it  is  best  to  consider  different  modifications  of  the 
standard  approach  separately  rather  than  in  combination. 

To  our  knowledge  there  is  only  one  convergence  result  in  the 
literature  for  multi-level  methods  on  locally  refined  grids.  This  is 
given  in  Bank  and  Dupont  (1978)  where  a  two-level  multi-level  iteration 
is  shown  to  converge  at  a  rate  independent  of  the  mesh  size.  Since  this 
result  relies  only  on  local  properties  of  the  finite  element  space  and 
not  on  global  approximation  properties,  it  applies  automatically  to 
locally  refined  grids.  The  limitation  to  this  result  is  that  it  deals 
only  with  their  two-level  scheme.  With  such  a  scheme,  as  one  decreases 


the  mesh  size,  the  dimensions  of  both  the  finer  and  coarser  finite 
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element  spaces  grow.  In  their  algorithm  the  linear  system  for  the 
coarser  finite  element  space  must  be  approximately  solved  once  per 
multi-level  iteration.  Whether  this  is  done  directly  or  by  one  of 
the  common  iterative  methods,  the  cost  of  approximately  solving  this 
linear  system  grows  faster  than  linearly  in  the  dimension  of  the 
linear  system.  As  a  consequence,  such  a  two-level  scheme  can  never  be  of 
optimal  order.  It  may,  however,  be  quite  attractive  as  a  practical 
algorithm  since  inversion  of  the  lower  dimensional  linear  system  for 
the  coarser  space  may  be  much  cheaper  than  the  cost  of  inverting  the 
linear  system  for  the  finer  space  in  which  the  solution  is  sought. 

The  use  of  direct  methods  to  solve  these  linear  systems  for  the  coarser 
space  is  particularly  attractive  since  after  the  first  multi-level 
iteration  only  back  solves  are  required. 

One  might  expect  that  a  simple  generalization  of  their 
convergence  proof  for  the  two-level  case  could  be  used  to  show  optimal 
order  convergence  of  a  multi-level  method  using  many  levels  on  locally 
refined  grids.  In  some  cases  it  can,  for  example,  for  Poisson's  equation 
on  certain  finite  element  grids.  In  general,  however,  the  error 
reduction  per  iteration,  y  £  (0,  1),  of  their  two-level  scheme,  while 
independent  of  the  mesh  size,  depends  strongly  on  the  PDE  and  on  certain 
properties  of  the  discretization.  The  inductive  argument  used  to  show 
the  optimal  order  convergence  of  multi-level  methods  using  many  levels 
requires  that  we  be  able  to  choose  a  sufficiently  small  value  of  y. 

For  example,  for  the  first  algorithm  in  this  chapter,  y  must  be  chosen 
in  (0,  -j-> .  Since  there  is  no  easy  way  of  controlling  y  in  their  two-level 
scheme,  it  seems  to  be  quite  hard  to  get  general  optimal  order  convergence 
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results  from  a  direct  generalization  of  their  two-level  result.  This  is 
unfortunate,  since  the  argument  showing  the  convergence  of  their 
two-level  algorithm  is  far  simpler  than  that  given  here. 

The  results  here  and  the  technique  of  their  proof  are  actually 

much  more  closely  related  to  the  second  kind  of  algorithm  considered  in 

their  paper.  This  second  algorithm  is  a  recursive  algorithm  using  many 

grid  levels  just  as  the  algorithm  in  the  last  chapter  did.  Their 

algorithm  is  shown  to  converge  in  0(N>  operations  on  quasi-uniform 

finite  element  spaces.  Unlike  the  results  in  chapter  three,  their 

results  apply  also  to  non-convex  domains,  on  which  the  elliptic  problem 
2 

is  less  than  H  regular.  The  analysis  of  algorithms  for  locally  refined 
grids  here,  which  also  applies  to  non-convex  domains,  is  largely 
modeled  on  their  analysis  of  the  second  algorithm  in  their  paper.  The 
notations  here  have  also  been  kept  almost  the  same  as  theirs. 

Analysis  of  multi-level  methods  is  harder  than  that  of  most 
iterative  methods  since  it  involves  both  algebraic  and  approximation 
questions.  The  algebraic  and  complexity  questions  for  locally  refined 
grids  are  slightly  harder,  since  the  dimensions  of  the  finite  element 
spaces  used  in  the  multi-level  algorithm  must  be  taken  into  account  more 
carefully.  By  contrast,  the  approximation  questions  for  locally  refined 
grids  seem  to  be  a  great  deal  harder.  In  the  first  half  of  this  chapter 
we  consider  the  algebraic  and  complexity  questions  involved  in  treating 
multi-level  methods  on  locally  refined  grids.  The  analysis  here 
largely  follows  that  in  the  paper  by  Bank  and  Dupont,  while  the  algorithms 
themselves  are  motivated  by  the  adaptive  finite  element  code  of  Bank 
and  Sherman.  Two  algorithms  are  considered;  the  first  is  similar  to 
that  in  the  Bank  and  Sherman  code.  The  second  is  a  modification  yielding 
an  improved  complexity  bound.  This  modification,  which  appears  to  be 


* 


« 


41 

new,  allows  greater  freedom  In  the  choice  of  finite  element  spaces 

used  and  may  have  significant  practical  applications. 

The  second  half  of  this  chapter  deals  with  the  approximation 

questions  involved  in  treating  locally  refined  grids.  The  approximation 

result  proved  there  is  the  main  result  of  this  thesis.  In  previous  work 

on  multi-level  methods,  the  approximation  questions  were  answered  by 

using  the  elliptic  regularity  of  the  problem  to  show  smoothness  of  the 

functions  being  approximated.  Then  the  approximation  error  could  be 

estimated  using  standard  finite  element  techniques  including  duality 

arguments  and  the  Bramble-Hilbert  lemma.  The  approximation  result  here 

uses  a  completely  different  kind  of  argument  making  no  use  of  the 

regularity  of  the  problem.  This  is  why  it  applies  to  non-convex 

2 

domains  where  the  regularity  is  weaker  than  H  .  Perhaps  of  greater 
interest  from  a  computational  point  of  view  is  the  fact  that  this  approx¬ 
imation  result  relies  only  on  local  properties  of  the  finite  element 
spaces  and  the  PDE  operator.  In  consequence, it  would  be  quite  easy  to 
compute  rigorous  upper  bounds  on  the  rate  of  convergence  of  the  algorithms 
in  this  chapter.  Up  until  now  there  were  good  heuristic  bounds  on  the 
rate  of  convergence  of  multi-level  methods  on  finite  difference  grids, 
Brandt  (1977),  and  in  a  few  cases  rigorous  bounds.  No  previous  bounds 
for  the  rate  of  convergence  of  multi-level  methods  on  irregular  finite 
element  grids  appear  to  be  available. 

4.2.  Notation 

To  develop  the  algorithms  considered  in  this  chapter,  it  is 
necessary  to  go  into  greater  detail  concerning  the  finite  element  spaces 
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used  than  was  required  in  the  last  chapter.  Additional  notation  must 
also  be  developed.  These  are  the  tasks  of  this  section.  To  keep  the 
mathematical  issues  in  this  chapter  clear,  our  analysis  will  be 
restricted  to  triangular  and  rectangular,  elements .  Though  the 
analysis  in  this  chapter  seems  to  be  extendable  to  all  or  nearly  all 
families  of  finite  elements,  the  details  of  the  proofs  here  can  be 
considerably  harder  for  some  families  of  elements.  Using  the  square 
and  triangular  elements  considered  here  permits  considerable  simplifi¬ 
cation  of  the  analysis  and  may  also  cover  the  cases  of  greatest  interest, 
since  the  majority  of  finite  element  computation  seems  to  be  done  with 
these  elements.  Linear  C°  triangular  elements,  called  Turner  triangles, 
are  used  in  the  Bank  and  Sherman  adaptive  mult'i-level  code,  so  at  least 
for  polygonal  domains,  the  solution  space  used  by  their  code  is  covered 
by  the  theory  here.  However,  as  will  be  seen, their  algorithm  violates 
the  assumptions  of  the  theory  developed  here  in  a  number  of  other  ways. 

There  are  really  two  reasons  why  more  detailed  specification 
of  the  finite  element  spaces  used  in  the  multi-level  algorithm  is  needed 
here  than  was  required  in  the  last  chapter.  The  first  is  that  locally 
refined  spaces  are  intrinsically  more  complex  than  quasi-uniform  spaces 
since  they  are  not  characterized  by  a  single  mesh  parameter.  The 
second  is  that,  unlike  the  analysis  in  the  last  chapter,  which  was 
based  on  standard  finite  element  error  estimates,  the  analysis  here  is 
based  on  an  approximation  result  developed  here.  The  proof  of  this 
result  depends  on  detailed  analysis  of  the  finite  element  spaces  used. 
While  this  result  could  have  been  proven  more  abstractly,  the  approach 
here  is  certainly  more  straightforward.  The  reader  with  a  background 
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in  finite  element  theory  should  have  little  trouble  in  seeing  the 


generalizations  to  elements  with  higher  order  continuity. 

The  two  elements  considered  here  are  the  standard  triangular 
and  square  elements.  The  C*“*  triangular  elements  of  degree  k,  k  _>  1, 
use  as  element  trial  space  the  space  of  polynomials  in  x  and  y 
of  degree  at  most  k.  The  square  elements  of  degree  k,  k  ^  1,  are 
tensor  product  elements  based  on  the  trial  space  0^  consisting  of  poly¬ 
nomials  whose  terms  are  of  degree  k  or  less  in  x  and  y  separately.  Both 
types  of  elements  may  be  taken  as  nodal  finite  elements  with  node 


placement  as  in  figures  1  and  2. 
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The  nodal  parameters  used  are  just  the  function  values  at  these  points, 
providing  a  substantial  simplification  of  the  analysis  in  this  chapter 
over  that  which  would  be  required  for  the  more  complex  Hermite  elements. 
Equally  important  from  our  point  of  view,  for  any  k  1 ,  either  type 
of  element  constitutes  an  affine  family,  Ciarlet  (1979).  What  this  means 
is  that  we  may  fix  a  reference  element  in  the  plane,  on  which  we  have 
a  reference  trial  space  L,  which  will  be  either  or  <2^  for  triangular 
or  square  elements,  respectively.  For  any  other  element  of  the  same 
type  anywhere  in  the  plane  we  can  find  an  affine  transformation  from 
the  reference  element  to  the  given  element.  This  affine  transformation 
carries  the  region  occupied  by  the  given  element,  the  nodes,  and  the 
element  trial  space  from  the  reference  element  to  any  other  element 
of  the  same  type  anywhere  in  the  plane. 

Affine  transformations  are  mappings  of  the  form 
x'  =  Ax  +  b 

where  b  is  a  fixed  coordinate  vector  and  A  is  a  nonsingular  2*2 
matrix.  For  triangular  elements  all  affine  transformations  may  be 
used,  carrying  the  reference  element  into  the  entire  six  parameter 
family  of  triangles  in  the  plane.  For  square  elements,  A  may  be  taken 
as  a  multiple  of  the  identity.  That  is,  there  is  no  need  for  rotation 
and  distortion  in  this  case. 

The  point  of  all  this  is  that  affine  transformations  and  the 
use  of  reference  elements  provide  an  elegant  simplification  of  the  theory. 
The  reference  element  and  its  element  trial  space  basically  contain  all 
properties  intrinsic  to  this  type  of  finite  element.  On  the  other 
hand,  the  affine  transformations  contain  all  information  about  the 
geometry,  mesh  size,  and  distortion  of  any  particular  element  in  the 


finite  element  grid.  Thus,  analysis  of  a  finite  element  space  is 
conveniently  factored  into  two  parts,  a  property  we  will  use  to 
advantage  in  section  4.5. 

Aside  from  the  choice  of  element,  two  types  of  grids  need  to 
be  considered.  These  may  be  conveniently  labeled  "aligned"  and 
"unaligned"  grids.  Aligned  grids  are  those  which  are  standard  in  finite 
element  theory,  where  each  edge  of  an  element  is  the  edge  of  another, 
or  part  of  the  boundary.  Unaligned  grids  are  those  where  the  edge  of  an 
element  may  be  only  part  of  the  edge  of  another.  See  figures  3,  4,  and  5 
for  examples  of  these  two  types  of  grids.  In  general,  triangular  grids 
may  always  be  chosen  to  be  aligned,  without  sacrificing  efficiency,  but 
it  is  often  necessary  to  choose  grids  of  square  or  rectangular  elements 
to  be  unaligned,  or  much  of  the  benefit  of  local  refinement  will  be  lost. 

To  be  admissible,  the  finite  element  spaces  must  be  C®  and 
satisfy  the  essential,  that  is  Dirichlet,  boundary  conditions.  On 
aligned  grids  there  is  no  problem  obtaining  continuity.  One  simply 
requires  agreement  of  the  nodal  parameter  on  the  edge  between  adjacent 
elements.  This  does  not  work  for  unaligned  grids.  To  simplify  this 
problem  we  restrict  the  unaligned  grids  as  follows.  Assume  that  every 
edge  of  an  element  is  either  part  of  the  boundary,  the  edge  of  another 
element,  or  the  union  of  the  edges  of  two  other  elements.  This  last 
situation  is  shown  in  figure  6. 

In  this  last  situation  we  consider  the  single  element  to 
dominate  the  two  elements  it  adjoins  across  this  edge.  For  convenience 
call  the  single  element  the  "larger"  since  it  almost  always  is  larger 
than  the  other  two.  Then  we  may  say  that  the  "larger"  dominates  the 


Figure  3.  Aligned  Triangular  Grid 


Figure  6.  Adjacent,  Unaligned  Elements 
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other  two  in  the  sense  that  its  nodes  on  this  edge  are  the  real  nodes 
of  the  finite  element  space  there.  The  nodes  of  the  "smaller"  elements 
may  be  called  parasitic  since  their  values  will  be  taken  from  the  larger 
element.  Parasitic  nodes  of  the  "smaller"  elements  in  figure  6  not 
coinciding  with  nodes  of  the  "larger"  element  are  labeled  with  a  "P." 

It  is  easy  to  see  that  this  system  preserves  the  required  continuity 
and  is  fairly  straightforward  to  code. 

One  more  item  needs  to  be  discussed  before  we  consider 
construction  of  the  locally  refined  spaces  needed  in  this  chapter.  This 
is  the  angle  condition.  It  was  mentioned  in  section  2.3  that  the  mass 
matrix  will  usually  be  well  conditioned  uniformly  in  the  level  of 
refinement  of  the  grid  as  long  as  the  element  geometry  does  not 
degenerate.  While  this  is  not  always  true  for  locally  refined  grids, 
it  is  still  important  that  the  element  geometry  not  degenerate.  One 
way  of  stating  the  required  condition  is  to  say  that  the  minimum  interior 
angle  of  any  triangle  in  the  grid  must  be  bounded  below  uniformly  in 

.  .  OO 

the  level  of  refinement.  More  formal lv,  if  (T.r.  ,  is  the  family  of 

J  J=1 

triangulations  being  considered,  we  want 

(4.2.1)  inf  {  min  <J>_}  >  0 

j>l  T£t. 

J 


where  <J>  is  the  smallest  interior  angle  of  the  triangle  T.  There  are 
several  alternate  formulations  of  this  condition  in  the  finite  element 
literature,  Oden  and  Reddy  (1976),  but  this  one  is  adequate  for  our 
purposes.  Obviously,  no  such  condition  is  needed  on  square  grids. 

We  next  turn  to  consideration  of  the  class  of  finite  element 
grids  on  which  the  multi-level  algorithms  of  this  chapter  can  be  shown 
to  converge  rapidly.  This  class,  which  includes  both  locally  refined  and 
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quasi-uniform  spaces,  is  quite  large  and  may  be  thought  of  as  an  indexed 
family  {M^}^  of  spaces.  It  consists  of  all  spaces  which  can  be 
generated  by  the  refinement  process  described  below  and  which  satisfy 
an  additional  condition  given  in  section  4.4.  An  0(N)  complexity  bound 
holds  uniformly  on  the  class  of  spaces  in  {M  }  .  satisfying  this 

condition  uniformly.  In  general,  this  class  of  spaces  with  a  uniform 
0(N)  complexity  bound  is  large  enough  to  show  the  improved  convergence 
to  singular  solutions  obtainable  by  discretization  with  locally  refined 
grids.  A  partial  exception  to  this  is  the  case  of  point  singularities, 
since  the  best  discretizations  for  these  problems  do  not  satisfy  the 
conditions  of  section  4.4  uniformly.  For  these  problems  0(N)  bounds 
can  be  shown  by  the  theory  here  only  for  suboptimal,  but  still  very  good 
locally  refined  grids. 

Fixing  a  coarse  space  M  ,  the  refinement  process  described 
below  will  generate  sequences  {M  ,  8  G  8  of  finite  element  spaces, 

P>J  3-1 

where  the  index  8  ranges  over  the  uncountable  family  of  choices  allowed 
in  this  refinement  process.  The  family  of  spaces  {M^}  ^  is  just  the 

union  of  all  spaces  which  can  be  produced  by  this  refinement  process, 
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The  multi-level  algorithms  here  actually  use  only  the  finite 
family  of  spaces  {AL  .}.  ,  for  fixed  3  e  B  and  fixed  k  >  2,  in  their 

3,  j  J=1  - 

computation.  Since  3  is  fixed,  it  is  convenient  to  drop  it,  although  we 

retain  j,  and  in  fact  consider  the  entire  sequence  of  refinements  {M  } 

3  J=1 

since  that  simplifies  the  statements  of  some  of  the  theorems  here.  One 
word  of  warning  is  needed,  however.  Though  the  spaces  in  the  sequence 

become  infinite  dimensional,  they  do  not  usually  become  dense  in 
Up  and  thus  convergence  to  the  true  solution  will  not  ordinarily  be 
obtained  in  {M  }  In  order  to  get  a  sequence  of  solutions  converging 

to  the  true  solution,  one  must  go  back  to  consideration  of  the  family 


{M 

u  o^A 


Description  of  the  families  {M  of  locally  refined  grids 


is  essentially  the  same  for  grids  based  on  square  and  triangular  elements. 

We  consider  first  the  triangular  case  since  it  is  the  case  of  greatest 

interest  in  this  chapter.  In  this  case,  the  family  {M.}.  ,  of  finite 

3  J=1 

element  spaces  will  be  based  on  a  nested  family  of  triangulations  {x 
By  nested  we  mean  that  each  triangle  T  of  T.,  j  _>  2 ,  is  contained  in  a 

,  .00 

triangle  of  Once  we  have  chosen  the  triangulations  and 

the  particular  element  to  be  used,  the  family  of  finite  element  spaces 

r  i  i  rx> 

lAl.i.  ,  is  completely  determined. 

.1  ,1  =  1 

,  ,  OO 

The  construction  of  the  nested  family  of  triangulations  IT  j.  - 

j  j=i 

may  be  done  by  a  process  Bank  and  Sherman  call  regular  subdivision.  In 
regular  subdivision  of  a  triangle  T,  it  is  divided  into  four  smaller  triangles 
by  pairwise  connecting  the  midpoints  of  the  edges  together  as  shown  in 
figure  7.  Let  be  a  fixed,  aligned  triangulation  of  fi.  Now  inductively 
suppose  T.  has  been  chosen  for  some  j  >  1.  Select  a  subset  t ^  of  .  We 


I 
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may  take  Tj+^  as  the  (unaligned)  triangulation  of  f2  formed  by  regular 
subdivision  of  the  triangles  in  while  leaving  the  remaining  triangles 
of  unaltered.  Proceeding  is  this  way,  the  entire  family  of 


triangulations  { T ^  } j _ may  be  constructed. 


Figure  7.  Regular  Subdivision  of  a  Triangular  Element 
In  our  approach  the  subsets  x  ^  C  x  ,  j  >.  1,  of  triangles 
to  be  refined  cannot  be  chosen  arbitrarily,  but  most  satisfy  certain 
conditions  necessary  for  the  operation  and  analysis  of  the  algorithms 
here.  Let  0  ^  c  f2  be  the  subdomain  convered  by  triangles  in  ,  for 

j  1.  Fox  simplicity  of  notation  it  is  convenient  to  set  =  f). 

The  first  condition  required  is 


(4.2.2) 


Vi  cVJiK 


This  condition  means  that  for  j  >  1,  only  triangles  of  Xj  which  were 
formed  by  regular  subdivision  of  triangles  in  T ^  ^  may  be  subdivided 
in  creating  Tj+j*  Equivalently,  for  j  >  1,  each  triangle  of  tj+j  Is 
contained  in  a  triangle  of  x^ . 

Now  assume  that  the  diameters  d^,  of  the  triangles  T  in  the 
coarsest  triangulation  X^  satisfy 


(4.2.3) 


-  hl  <  dT  <  o  hl  , 


for  some  fixed  a  >  1  and  mesh  parameter  h^.  As  before  set 


4 
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(4.2.4)  h.  =  P1_j  h1  ,  j  >  1  , 

where  p  is  a  mesh  ratio,  which  for  the  regular  subdivision  process  here 
is  two.  Conditions  (4. 2. 2)- (4. 2 . 4)  imply  that  the  diameters  of 
triangles  in  T ^ ,  j  _>  1,  must  satisfy 

(4.2.5a)  -h.  <d  <oh.  forTC  jj. 

O  j  -  T  -  j  j 

(4.2.5b)  —  h.  <  d_  <  u)  h .  for  T  C  j).  \ 

ax—  T—  i  i  i+I 

with  1  <_  i  <_  j  -  1  . 

The  overbar  here  denotes  closure,  necessary  since  triangles  are  closed, 

r  .oo 

and  our  domains  1S2.}.  _  are  not. 

j  j=i 

Intuitively,  the  triangles  in  are  all  of  about  the  same 
size.  By  (4.2.3)  or  (4.2.5a)  their  diameters  are  all  comparable  to  h^. 
The  triangles  in  come  in  two  sizes:  fine  ones  in  with  diameter 
comparable  to  h^,  produced  by  regular  subdivisions,  and  coarse  ones 
outside  (2 2  with  diameters  about  h^.  In  general,  may  be  thought  of 
as  having  j  sizes  of  triangles,  the  finest  lying  in  12^.  Of  course, 
this  intuitive  point  of  view  should  not  be  pushed  too  far.  Because  of 
the  constant  0  in  (4.2.5)  these  sizes  may  overlap  a  great  deal. 

Only  one  other  condition  is  really  essential  for  triangular 
grids.  That  is,  that  the  number  of  elements  adjoining  a  given  element 
along  one  edge  be  at  most  two.  The  "two"  here  is  arbitrary,  and  any 
other  finite  bound  could  be  chosen,  but  in  practice  two  is  probably  all 
one  wants.  The  coding  is  difficult  enough  as  it  is.  This  restriction, 
that  at  most  two  "smaller"  elements  adjoin  one  "larger"  element  across 
an  edge,  together  with  the  angle  condition  (4.2.1)  limits  the  rate  of 
change  of  the  mesh  size  as  one  can  easily  show. 
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Locally  refined  grids  of  square  elements  may  be  constructed  along 
the  same  lines  as  the  triangular  grids  just  considered.  One  begins  with 

an  aligned  square  grid  G^  on  ft  and  then  generates  a  family  of  finer 

0° 

(unaligned)  grids  {G  }j_j,  using  an  analogous  regular  subdivision  process. 

f  "I  00 

This  process  is  shown  in  figure  8.  The  subdomains  {ft,}.  ,  can  be  defined 

j  j=i 

analogously,  and  the  analogs  of  (4 . 2. 3)- (4 . 2 . 5)  will  hold,  now  with  the 
notation  "d^"  instead  of  "d^",  using  "e"  for  "element."  Since  such 
square  grids  are  only  being  considered  in  this  section  and  will  thereafter 
be  dropped  for  simplicity,  there  is  no  need  to  go  into  these  details 
more  deeply. 


Figure  8.  Regular  Subdivision  of  a  Square  Element 
A  problem  arises  on  locally  refined  grids  not  present  on 
quasi-uniform  grids.  Unless  special  care  is  taken,  the  condition  numbers 
of  the  mass  matrices  for  a  family  of  locally  refined  finite  element  spaces 
may  not  be  uniformly  bounded.  For  aligned  quasi-uniform  grids  and  the 
usual  types  of  elements,  including  those  considered  here,  the  uniform 
boundedness  of  the  condition  number  of  the  mass  matrices  is  equivalent 
to  the  angle  condition  (4.2.1).  Unaligned  quasi-uniform  grids  are  rarely 
considered,  for  obvious  reasons,  but  the  same  thing  can  be  shown  for 


such  grids. 


The  problem  with  locally  refined  grids  is  that  the  number  of 
basis  functions  which  overlap  with  a  given  basis  function  may  be  very 
large,  becoming  unbounded  as  the  grid  is  refined.  This  violates  the 
hypotheses  of  the  result  of  Fried  described  in  section  2.3,  and  in  fact 
mass  matrices  with  condition  numbers  becoming  unbounded  do  occur.  A 
simple  finite  element  grid  where  this  problem  occurs  is  shown  in 
figure  9.  Here  bilinear  elements  are  used  with  the  standard  "hat  function" 
or  "pagoda"  basis  functions.  Notice  that  the  conditioning  problem  in 
this  example  occurs  despite  the  fact  that  the  mesh  size  changes  slowly 
in  the  sense  that  the  diameters  of  adjacent  elements  differ  at  most  by 
a  factor  of  two.  Proper  scaling  of  the  basis  functions  mitigates  the 
problem  but  does  not  eliminate  it. 


±± 

Figure  9.  Finite  Element  Space  with  Poorly  Conditioned  Mass  Matrix 

The  only  real  solution  is  to  arrange  things  so  that  each  basis 
function  overlaps  only  a  bounded  number  of  elements,  with  the  bound 
independent  of  the  level  of  refinement.  This  can  be  done  either  by 
modifying  the  basis  elements  to  prevent  excessive  overlap,  or  by  selecting 
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only  grids  where  such  overlap  does  not  occur.  The  second  approach  is 
probably  preferable  in  most  cases  since  it  tends  to  be  simpler.  Consider 
the  approach  of  modifying  the  basis  functions  for  the  simple  bilinear 
elements.  The  situation  is  similar  for  the  other  square  and  triangular 
elements  here.  Ordinarily,  the  nodal  parameter  at  the  intersection  of  the 
four  elements  shown  in  figure  10  would  be  associated  with  the  "hat 
function"  whose  support  is  the  union  of  the  four  elements.  When  one  or 
more  of  these  four  elements  is  refined,  the  basis  function  associated  with 
its  center  node  can  be  modified  to  have  support  on  a  smaller  region. 

Up  to  rotation,  there  are  four  cases  to  consider,  as  shown  in  Figure  11, 
Clearly  the  case  where  all  four  elements  are  refined  need  not  be  dealt 
with  since  then  the  same  considerations  apply  on  the  resulting  smaller 
elements.  The  dark  lines  in  figure  11  show  the  boundaries  of  the 
supports  of  the  modified  basis  functions.  The  squares  marked  with 
asterisks  may  be  subject  to  further  refinement  without  violating  the 
requirement  that  the  mesh  size  change  slowly.  The  point  is  that  all 
squares  so  marked  lie  outside  the  support  of  these  basis  functions,  so 
the  number  of  elements  these  basis  functions  overlap,  and  hence  the 
number  of  other  basis  functions  they  overlap,  will  be  bounded  by  a  fixed 
constant  independent  of  the  level  of  refinement.  The  only  difficulty  with 
this  approach  is  that  computations  of  the  element  integrals  is  somewhat 
more  complicated  since  these  basis  functions  are  more  complex.  This 
approach  was  suggested  to  the  author  by  D.  Gannon,  who  has  done  considerable 
work  with  rectangular  tensor  product  elements  on  locally  refined  grids. 


Gannon  (1980). 
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Figure  LI.  Modified  Basis  Functions 
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In  the  second  approach,  one  retains  the  simple  Lagrangian  basis 
functions,  "hat  functions,"  and  their  higher  order  analogues,  but 
restricts  or  modifies  the  grids  so  that  the  problem  of  one  basis  function 
overlapping  an  increasingly  large  number  of  others  never  arises.  With 
triangular  grids  modification  is  not  difficult.  This  is  the  approach 
used  in  the  Bank  and  Sherman  code.  Beginning  with  an  unaligned 
triangulation  T.,  1  <_  j  <  00 ,  belonging  to  a  nested  family  of  triangulations 
satisfying  the  conditions  given  so  far,  their  algorithm  constructs 
an  aligned  triangulation  ?  ,  refining  T  a  on  which  this  condition  number 
problem  does  not  arise.  This  is  done  in  two  steps.  In  the  first  step 
each  triangle  with  the  property  that  there  are  vertices  of  other 
triangles  in  the  middle  of  more  than  one  of  its  edges  is  regularly 
refined,  as  in  figure  12.  Once  all  such  triangles  are  refined,  including 
those  produced  in  carrying  out  this  process,  the  only  triangles  in  the 
grid  are  those  with  either  one  vertex  in  the  middle  of  an  edge  or  no 
vertices  in  the  middle  of  its  edges.  The  former  are  refined  by  a 
different  subdivision  process,  called  "green"  refinement,  shown  in 
figure  13.  The  resulting  aligned  triangulation  T  is  then  used  to  generate 
the  finite  element  space  used  by  the  multi-level  algorithm. 

Notice  that  green  refinement  may  decrease  the  smallest  interior 
angle  present  on  the  grid.  This  would  be  a  problem  if  successive  green 
refinements  were  applied.  In  Bank  and  Sherman's  approach  this  is  avoided 
by  first  generating  the  unaligned  triangulations  {t  }  ^  and  then 

forming  the  corresponding  aligned  triangulation  rather  than 

taking  Tj+1>  j  >  1.  as  a  refinement  of  t  .  This  ensures  that  the  angle 
condition  (4.2.1)  is  satisfied.  Even  so,  the  triangulations  it  }  ^  used 

by  their  code  are  not  generally  admissible  in  the  theory  here,  even  when 
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Figure  12.  Regular  Refinement  of  Triangles  with  Edges  Containing  Vertices 


Figure  13.  Green  Refinement  of  a  Triangle 
the  underlying  unaligned  triangulations  {t  }  ^  satisfy  all  of  our 

assumptions.  The  problem  is  that  in  producing  an  aligned  triangulation 
.  j  >  1,  a  triangle  may  be  green  refined  while  in  a  finer  unaligned 
triangulation  T\,  i  >  j,  this  same  triangle  may  be  regularly  refined. 

If  so,  the  triangulations  1 T ^  / j _ will  not  be  nested,  nor  will  the  resulting 
finite  element  spaces  {M^ }  ^.  While  this  does  not  seem  to  be  a  problem 

for  their  code,  it  is  certainly  a  problem  for  the  theory  here.  The  cure 
for  this  seems  to  be  to  require  the  unaligned  triangulations  {t^},  ^ 
to  contain  only  triangles  with  a  vertex  of  another  triangle  in  at  most 
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r  t  00 

one  of  their  edges.  If  so,  in  going  from  l t ^  i j _ to  the  corresponding 
aligned  triangulations  only  green  refinement  will  be  done.  Since 

the  triangles  of  ,  j  ^  1,  green  refined  in  producing  T  must  lie 
outside  fij  (the  portion  of  lying  in  is  aligned),  there  is  no 
danger  of  subsequently  regularly  refining  them.  Further,  it  is  easy 
to  see  that  those  triangles  of  f  produced  by  green  refinement  must  be 
green  refined  in  producing  T  ,  i  >  j,  so  the  finite  element  spaces  based 
on  the  family  of  triangulations  It)  ^  will  be  nested  and  the  theory  of 
this  chapter  will  go  through  for  them. 

The  advantage  of  this  approach  is  that  the  condition  number 
of  the  mass  matrix  for  a  finite  element  space  based  on  an  aligned 
triangulation  will  be  uniformly  bounded,  as  long  as  the  angle  condition 
is  satisfied.  This  is  so  for  linear  elements  since  each  pyramid  basis 
function  overlaps  only  those  elements  meeting  at  the  vertex  with  which 
its  nodal  parameter  is  associated.  The  support  is  thus  a  polygon  formed 
from  these  elements,  as  in  figure  14.  The  other  kinds  of  basis  functions 
required  for  higher  order  elements,  those  associated  with  nodes  on  the 
edges  or  interiors  of  triangles,  are  not  a  problem  either,  since  their 
supports  are  smaller  than  those  associated  with  the  vertices. 


Figure  14.  Support  of  a  Basis  Function  on  an  Aligned  Triangu1 ation 
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Though  the  above  approach  is  quite  attractive  for  triangular 
grids,  since  aligned  triangulations  are  convenient  in  programming,  it 
does  not  generalize  directly  to  square  elements,  since  there  is  no 
analog  of  green  refinement  for  these  elements.  Instead,  we  must  either 
use  the  modified  basis  functions  described  above,  or  restrict  the  rate 
of  change  of  the  mesh  size  so  that  the  conditioning  problem  does  not 
occur  when  using  the  standard  basis  functions.  The  simplest  way  to 
describe  the  required  restriction  on  the  rate  of  change  of  the  mesh 

,  .00 

size  is  in  terms  of  the  process  for  constructing  the  grids  of 

square  elements  we  wish  to  consider.  Say  that  two  elements,  e^,  £  G., 

1  £  j  <  00 ,  touch  if  they  have  any  point  in  common;  that  is,  if  they 
share  a  vertex  or  part  of  an  edge  as  in  figure  15.  Let  Gj  ,  j  _>  1,  be 
the  set  {e  €  Gj  :  e  C  That  is,  Gj  is  the  set  consisting  of  the 

most  refined  elements  of  G  ^ .  Define  the  interior  of  Gj ,  denoted 
int(Gj)  as  the  set  of  elements  of  Gj  not  touching  any  elements  in 
Gj \  Gj .  Then  a  sufficient  condition  for  the  mass  matrices  for  the  locally 

r  i  t 00 

refined  spaces  based  on  the  grids  {G^  > ^ to  be  uniformly 

well  conditioned  is 

G^  C  int(Gj )  ,  j  >  1  . 

This  condition  suffices  for  the  mass  matrices  to  be  uniformly  bounded 


since  no  basis  functions  centered  at  nodes  outside  of  Gj  have  supports 
extending  into  int(Gj).  It  follows  that  they  only  overlap  elements  outside 
SI  ,  and  so  they  only  overlap  a  bounded  number  of  elements. 
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Figure  15.  Touching  Elements 
With  the  two  types  of  elements  considered  here  and  the 
variety  of  ways  of  ensuring  the  uniform  boundedness  of  the  condition 
numbers  of  the  mass  matrices,  a  number  of  options  are  available. 
However,  it  would  be  unwieldy  to  treat  the  theory  in  the  rest  of  this 
chapter  in  the  same  generality,  particularly  the  work  in  section  4.5, 
which  relies  on  detailed  analysis  of  the  finite  element  spaces  used. 
For  this  reason,  we  restrict  the  analysis  that  follows  to  C®  spaces 
on  aligned  triangular  grids.  This  choice  is  motivated  partly  by 
convenience,  and  partly  by  the  use  of  such  spaces  in  the  Bank  and 
Sherman  code.  There  would  be  no  difficulty  in  treating  the  other 
types  of  finite  element  spact  ;  considered  here,  however,  the  interpo¬ 
lation  mappings  of  section  4.5  are  slightly  more  complicated  to 
analyze  on  unaligned  grids. 
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4.3.  Algorithms 

In  this  section  we  describe  the  multi-level  algorithms  to  be 
considered  in  this  chapter  and  begin  the  analysis  of  their  convergence. 

The  first  of  these  algorithms  is  identical  to  one  of  those  described  by 
Bank  and  Dupont  for  quasi-uniform  grids  and  is  similar  to  algorithms 
described  by  Nicolaides,  Hackbusch,  Brandt,  Bokhvalov,  and  others.  It 
is  also  similar  to  the  algorithm  used  in  the  adaptive  finite  element 
code  of  Bank  and  Sherman  though  there  is  a  significant  difference  here 
which  will  be  discussed.  The  second  algorithm  to  be  considered  here 
appears  to  be  new.  It  can  be  used  to  show  the  0(N)  complexity  bound 
given  in  the  next  section  under  weaker  assumptions  than  otherwise  necessary. 

Both  of  the  algorithms  described  here  differ  substantially 
from  that  described  in  the  last  chapter.  There,  to  approximately  solve 


the  finite  element  equations  on  a  member  M.,  j  >  2,  of  a  family  {M.} 


j 
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of  finite  element  spaces,  approximate  solution  of  a  related  problem  in 
M.  -j  was  required  twice;  once  at  the  beginning  to  obtain  a  starting 
value  for  smoothing  iteration  on  M. ,  once  at  the  end  of  the  iteration  as 
a  final  correction.  The  algorithms  in  this  chapter  and  other  multi-level 
algorithms  in  the  literature  use  approximate  solution  of  the  equations 
for  M.  j  more  often.  After  obtaining  a  starting  value  from  Mj_^,  the 
accepted  solution  is  obtained  through  a  sequence  of  multi-level  (outer) 
iterations.  F.ach  of  these  outer  iterations  consists  of  a  number  of 
smoothing  (inner)  iterations  on  M  ,  followed  by  approximate  solution  of 
a  related  problem  in  Mj  Thus,  approximate  solution  of  the  related 

problem  in  M  ^  is  required  an  indeterminate  number  of  times  depending 
on  how  far  the  outer  iteration  on  is  carried. 
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This  kind  of  algorithm  is  somewhat  more  complicated  than  that 

of  chapter  three.  However  it  made  sense  to  analyze  the  application  of 

this  kind  of  algorithm  to  locally  refined  grids  rather  than  the 

application  of  algorithms  like  those  in  the  last  chapter  for  several 

reasons.  The  main  one  was  the  desire  to  establish  the  optimal  order 

convergence  of  an  algorithm  as  close  as  possible  to  those  being  used 

in  practice.  This  was  more  fully  accomplished  than  it  would  have  been 

if  we  had  examined  algorithms  more  like  those  in  chapter  three.  It 

might  also  be  pointed  out  that  relatively  little  of  the  complexity 

of  this  chapter  is  caused  by  the  choice  of  algorithm. 

Description  of  the  algorithms  here  and  analysis  of  their 

convergence  properties  will  involve  the  eigenfunctions  and  eigenvalues. 

Following  the  notation  in  chapter  three  let  Nj  be  the  dimensions  of  , 

1  £  .1  <  «°.  Then  for  the  problem  (2.1.3)  we  have  eigenfunctions 

in  M.  and  eigenvalues  satisfying 

l  i=l  j  i  i=l  7  b 

<t>)  =  <f>),  4>  e  M  . 

Since  more  than  one  finite  element  space  may  be  involved  in  our 

considerations,  we  will  retain  the  superscript  j.  As  in  chapter  three 

2 

we  may  replace  the  L  inner  product  by  the  bilinear  form  b.(*,  •)  defined 
there,  where  we  now  use  the  subscript  j  to  keep  track  of  the  finite 

.  QO 

element  space  in  (M^ }  ^  giving  rise  to  this  bilinear  form.  With  this 

replacement  we  have  as  before  eigenfunctions  and  eigenvalues 

^i'^i-1  satl-sfyln8 

♦)  “  Xjj)  b  d>),  <t>  G  M  . 

As  usual  we  assume  the  eigenvalues  are  numbered  in  order  of  increasing 
magnitude: 


« 


f 


o  <  <  x^  <...<  )  , 

I  ~  l  —  Nj 

0  <  x(j)  <  x^j)  <...<  x£J)  . 

J 

For  both  of  the  algorithms  considered  in  this  chapter  the  basic 
multi-level  iteration  is  the  same.  In  both  cases  each  outer  iteration 
consists  of  a  sequence  of  smoothing  inner  iterations  followed  by 
approximate  solution  of  a  related  problem  in  a  coarser  finite  element 
space.  The  difference  between  the  two  algorithms  will  lie  in  the  manner 
of  solving  these  related  problems.  In  order  to  obtain  a  solution  in  a 
member  M^,  k  >  2,  of  a  family  of  locally  refined  spaces, 

multi-level  iteration  will  need  to  be  applied  to  all  of  the  spaces 
{ Afj  } j _ 2 »  thus,  the  iteration  will  be  described  for  any  space  in 

(Vl=r 

Beginning  with  a  trial  solution  u„ ^  in  M.  the  outer  iteration 

O.n  j 

produces  a  sequence  of  iterates  ^ug^n^n-l  conver8:tng  to  the  true 

discrete  solution  u  ^  ^  in  .  Each  of  these  outer  iterations  is 

similar  to  the  algorithm  of  the  last  chapter.  Beginning  with  a  trial 

solution  u<j)  G  M,  an  outer  iteration  produces  an  improved  solution 
0,n  j 

(j) 

U0,n+1  given  by 


u(3)  =  u(1) 

0 ,n+l  m+l,n  ’ 

where  {uph!>+|  is  a  sequence  of  inner  iterates.  All  but  the  last  of 
?.,n  Jt=l 

these  inner  iterates  are  generated  by  a  simultaneous  displacement 

smoothing  iteration  on  M  .  The  last  uses  a  correction  based  on  the 

approximate  solution  of  a  related  problem  in  the  coarser  finite  element 

space  M  .  More  precisely,  beginning  with  a  trial  solution  u^^  6  M. , 
1  - 1  0 ,  n  j 

an  outer  iteration  consists  of  the  steps: 
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1.  For  £  =  1,...,  m  define  u.  £  M.  by 

£,n  2 

(4.3.1)  -  “juj.n.  *>  “ 


ryjy  (afu^  n»  4>)  -  (f.  40)  »  4>  e  M 


Letting  v  £  M.  ^  be  given  by 


(4.3.2) 


a(v,  <p)  =  a(u  ,  40  -  (f,  <p)  ,<peM 

lu ,  n  j  “i- 


let  v  be  an  approximation  to  v  and  set 


(4.3.3) 


u  . .  =  u  -  V 

m+1 ,  n  m ,  n 


The  quality  required  in  the  approximation  of  v  by  v  will  be  specified 
in  theorems  4.1  and  4.3. 

The  intuitive  meaning  of  these  two  steps  is  that  in  the  first 
the  more  oscillatory  components  of  the  convergence  error  are  rapidly 
reduced  by  the  smoothing  iteration.  In  the  second  step  the  smoother 
components  of  the  error,  for  which  the  simultaneous  displacement 
iteration  is  ineffective,  are  reduced  through  approximate  solution  of  the 
related  problem  (4.3.2).  Consequently  convergence  will  be  rapid  for  all 
components  of  the  error. 

To  make  these  considerations  precise,  let  9  £  (0,  1)  be  an 

arbitrary  parameter.  Let  J(6j)  be  the  integer  such  that 

A.(J)  <  0,  A<j)  for  i  <  J(0  )  , 

1  “**  1  N .  —  1 

J 

x(j)  >  e  xf,j)  for  i  >  J(9  )  . 

1  j  Nj  j 


Similarly  let  J(6^)  be  the  integer  such  that 


« 


f 


xfJ)  <  e.  for  i  <  J(0.)  , 


i  -  j  N 


L(j)  >  0.  xf,j)  for  i  >  JT  (6  )  . 
i  j  N.  j 

3 

Then  we  can  define  the  subspace  S .  0 .  3.  0.  of  M.  by 

J  J  J  3  3 


=  span  {i|>^  ,  i  <  J(0  J}  , 


0.  =  span  i  >  j(0.)}  , 


5.  =  span  i  <  j(o . > }  , 


0^  =  span  {tfr-1  ,  i  >  J(8J}  . 

The  spaces  S.  and  S.  consist  of  relatively  smooth  functions  while  0. 

3  3  3 

and  0^  contain  more  oscillatory  functions.  For  any  0  €  (0,  1)  one  has: 

"j  ' s:  ©  ‘  ©  sj 

1 1.  fact,  0.  is  the  orthogonal  complement  of  in  M  in  both  the  energy 
and  L*-  inner  products,  while  ()_.  is  the  orthogonal  complement  of  in 
energy  and  the  inner  product  b^  (•,  •). 

For  the  time  being  only  the  subspaces  and  0^  will  be  important 
in  analyzing  the  convergence  of  the  multi-level  iteration  (4. 3.1)— (4. 3. 5) , 
since  the  approximation  question  involved  will  be  deferred  until  sections 
4.5-4. 7.  There  and  0^  will  play  an  important  role.  To  establish  the 
rapid  convergence  of  iteration  (4 . 3 . 1 )- (4 . 3. 3)  we  will  show  that  error 
components  in  the  oscillatory  space  0^  are  rapidly  reduced  in  step  1, 
while  components  in  the  smooth  space  are  reduced  in  step  2. 

The  multi-level  iteration  (4. 3. l)-(4. 3. 3)  differs  from  that 


in  the  Bank  and  Sherman  adaptive  finite  element  code  primarily  in  the 


smoothing  inner  iteration  used.  Their  code  uses  a  symmetric  Gauss-Seidel 


* 
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smoothing  iteration  rather  than  the  simultaneous  displacement  iteration 

used  here.  While  comparisons  for  quasi-uniform  grids  have  shown  that 

the  choice  of  smoothing  iteration  is  not  very  important,  Brandt  (1977), 

Bank  and  Sherman  (1979),  for  locally  refined  grids  the  choice  is  significant. 

In  Fourier  terms,  the  iteration  (A. 3.1)  rapidly  annihilates  Fourier  error 

components  whose  wave  length  is  comparable  to  h^ ,  the  finest  mesh  spacing 

of  M..  Away  from  ,  the  most  refined  region,  this  iteration  has  little 

effect.  The  reason  for  this  is  that  the  residual  on  the  right  hand 

side  of  4.3.1  is  multiplied  by  the  constant  factor  -  ttTTT  .  By  contrast 

\  (J  > 

Nj 

in  Jacobi  or  Gauss-Seidel  iteration  each  component  of  the  residual  is 

multiplied  by  the  inverse  of  the  corresponding  diagonal  entry  of  the 

stiffness  matrix.  In  other  words,  in  these  latter  cases  the  scaling 

2  2 

factor  looks  like  -ch,  .  rather  than  -ch . ,  where  h,  .  is  a  function 

local  j  local 

which  takes  on  a  different  value  for  each  component  of  the  residual 

vector,  depending  on  the  size  of  the  corresponding  elements.  As  a 

result  of  this  difference,  the  smoothing  iteration  used  in  the  Bank  and 

Sherman  code  is  effective  everywhere  on  the  grid,  rather  than  just  on  the 

most  refined  parts.  This  is  important  in  their  code  since  they  a ' low 

any  element  of  M  j  >  2,  to  be  refined  in  creating  a  finer  space  M  . 

J- 1 »  J 

We  do  not,  and  the  simultaneous  displacement  iteration  considered  here 
is  perfectly  adequate  in  this  simpler  case. 

The  choice  of  smoothing  iteration  is  intimately  connected  with 
the  decomposition  of  M  into  smooth  and  oscillatory  components.  The 
decomposition  here,  ©  ,  is  natural  only  for  simultaneous 

displacement  iteration.  The  hard  part  of  the  analysis  here  is  the 
approximation  question  involved  in  demonstrating  the  effectiveness  of 


« 


fi0-. 


9 


step  2  of  the  multi-level  iteration  (4 . 3. 1)- (A. 3 . 3) .  For  the  decomposition 

M,  =  S,  0.  and  the  case  considered  here  where  only  the  finest  elements 
j  j  v~/  J 

of  M  ^  are  refined  in  creating  ,  this  approximation  qu  stion  can  be 
answered.  The  author  was  unable  to  handle  the  more  general  case  where 
any  element  of  ^  may  be  refined. 

Though  this  approximation  result  will  not  be  shown  until 
section  4.6,  we  now  state  it  since  it  is  needed  for  the  analysis  in 
this  section.  Let  be  a  family  of  locally  refined  spaces 

satisfying  the  assumptions  of  section  4.2. 

Let  j  _>  2  and  e  (0,  1)  be  arbitrary,  and  for  any 
n  6  Sj  let  r)  £  ^  be  given  by 

a(n  -  n,  40  *  o  ,  <t>  e  M  , 

that  is,  h  is  the  elliptic  projection  of  D  onto  M ^ ^ .  Then 


(4.3.4) 


|n-n|||  <c0  ej'*  lllnlll  , 


where  c^  >  0  is  independent  of  j,  0^ ,  and  r).  This  inequality,  which 

is  the  main  result  of  this  thesis,  shows  that  eigenfunctions  in 

Mj  whose  eigenvalues  are  relatively  small  can  be  well  approximated  in  the 

coarser  finite  element  space  M. 

j-1 

With  these  preliminaries  behind  we  are  ready  to  consider  the 
first  convergence  theorem  for  the  multi-level  iteration  (4. 3. l)-(4. 3. 3) . 
This  theorem  is  the  same  as  one  found  in  Bank  and  Dupont  for  the  case 
of  quasi-uniform  grids  and  is  similar  to  theorems  proved  by  Nicolaides, 
Hackbusch,  and  others.  Our  notation  largely  follows  Bank  and  Dupont 
although  there  are  some  differences. 


i 
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Theorem  4.1.  There  exists  y  e  (0,  1)  and  a  positive  integer  m,  both 
independent  of  j  ^  2,  such  that  if  the  approximate  solution  v  of 
(4.3.2)  satisfies 


(4.3.o) 

then 

(4.3.6) 


v  -  v  <  y 


lUnH-l,n  “  1  Y  IK.n  "  U 


<J)| 


Except  for  our  use  of  (4.3.4)  rather  than  a  similar  inequality  for 
quasi-uniform  grids,  the  proof  of  this  theorem  is  identical  to  the 
proof  of  the  analogous  theorem  for  quasi-uniform  grids  given  by  Bank 
and  Dupont. 

Proof.  Let  0^  €  (0,  1)  be  free  to  be  chosen  later  and  consider  the 
decomposition  of  : 


Let  ,  0  _<  A  <_  m  +  1,  be  the  convergence  error. 


e0  =  u„  -  u 
Z  Z  ,n 


(j) 


since  the  second  subscript  on  u  plays  no  role  in  the  theorem.  Then 

GZ  =  n£  +  ’ 

for  some  e  5  and  £  e  0.  Expanding  e^  in  eigenfunctions, 


T.  a 
i=l 


.(j) 


i  i 


we  have  for  Z  m  0,  1,...,  m. 


:(J) 


N  Ai  l  ~ (1 ) 

£  V1  -  ^ 

11  AN 


« 


by  (4.  3. 1)  .  Thus  , 


Ills  III2  =  "  ai(1  "  ^7jT)2m  !H^j)ill2 

m  i=J(0)+l  1  JlV'  1 

N 


<  (1  -  0  )2m  E  a2  |j|^j)]||2 
J  i=J (0)+l 


iiic.ui  i  a  -  ym  Hicoiii  i  a  -  ym  ill e0 in  . 

IllnJH  £  |||n0l||  l  II! e0 III  • 

Let  X  and  f|  be  the  elliptic  projections  on  M.  ,  of  £  and  n  >  respectively 

j-i  mm' 

a(l  -£.♦)-  0  ,♦  e  M  , 

m  J-l 

a  (n  -  n.  0)  *  o  ,  <p  e  M  . 

m  j-l 

Since  X  is  the  projection  relative  to  the  energy  inner  product 


ms  -  yii  i  myii  i  a  -  o.)m  in eQ i 


and  by  (4.3.4) 


IS  -  njll  1  c08»4  UlnJII  i  c0e^  |||.0| 


Since  v  =  n  +  f  and  e  =  n  +  £  , 

m  m  m 


l«.  -  sill  ico8j1/4  III e0 III  +  (1  -  e.>*  nielli  . 


Then  since  e  =  e  -  v, 
nrfl  m  9 


|enrt.1lll  £  !!!=„  -  v|||  +  |||v  -  v| 


£  c„  5/4  |||.0IH  -8  (1  -  V-  III e0 III  ty2  lllvlll 

by  the  hypothesis  of  the  theorem.  Since  v  is  the  projection  of  em  on 


« 


t 
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e 

1  m1 


'0' 


Thus , 

«.3.7)  ||lVll 

To  complete  the  proof ,  let  y  £  (0,  1)  be  such  that 

2 
Y 


<  (c^174  +  (1  -  Sj)™  +  y2) 


2  <1. 


Then  choose  0.  £  (0,  1)  such  that 
1 


and  finally  m  such  that 


c.0 .  <* 
0  1  3 


(1  -  0  )m  <  |  . 


Since  these  choices  of  y>  6^  and  m  are  Independent  of  j,  the  proof 
is  complete.  □ 

Inequality  (4.3.5)  in  this  theorem  is  essentially  an  induction 
hypothesis.  Suppose  that  for  fixed  Jq  2  using  iteration  (4 . 3. 1)- (4 . 3. 3) 
with  j  =  j  causes  the  energy  norm  of  the  convergence  error  to  be  reduced 
by  the  factor  y  e  (0,  1)  according  to  this  theorem.  Also  suppose  we 
have  an  arbitrary  iteration  on  M.  .  which  reduces  the  energy  norm  of  the 

Jo_i 

convergence  error  on  M,  _.  by  the  same  factor  y  £  (0,  1)  at  each  iteration. 

J0 

Then  this  iteration  on  M  can  be  used  to  approximately  solve  the 

j0~ 1 

related  problem  (4.3.2)  on  M  .  arising  in  using  the  multi-level 

J0 

iteration  (4. 3. 1 )— (4 . 3. 3)  on  M  .  Beginning  with  the  trial  solution 

J0 

vo  -  0 

/  1°° 

in  M.  ,  the  iteration  on  M  produces  a  sequence  of  iterates  IvJ.  n 
J0-l  Jq-1  16  ^ 

satisfying 


•  v£  -  vl 


4 


9 


Thus,  two  of  these  iterations  on  M.  _  would  be  enough  to  satisfy  the 

J0_1 

hypothesis  (4.3.4)  of  theorem  4.1. 

According  to  this  theorem  the  constant  y  G  (0,  1)  by  which  the 
energy  norm  of  the  convergence  error  is  reduced  per  iteration  of 

(4 . 3 . 1) - (4 . 3.  3)  is  independent  of  j  _>  2.  It  follows  that  if  j^-l  2 
the  iteration  on  M,  ,  can  be  taken  as  the  multi-level  iteration 

Jo"1 

(4. 3.1)  —  (4. 3. 3)  with  j  *  j^-l.  ■*'n  case  we  must  solve  in  turn 

related  problems  on  M.  Again  this  can  be  done  recursively  with 

J0" 

iteration  (4. 3. l)-(4. 3. 3)  if  j  -2  >  2. 

To  make  these  ideas  clear,  we  may  consider  the  following 
informal  procedure.  Calling  this  recursive  procedure  with  j  = 
causes  n  of  the  multi-level  iterations  (4.3.1)  —  (4.3.3)  to  be  carried 
out  on  H.  . 


Procedure  Outer_l_on_M.  (n) 

I f  (j  =  1)  {solve  the  equations  on  directly;  then  return } 
else  (repeat  n  times{ 

1.  do  m  smoothing  inner  iterations  on  M  according  to  (4.3.1). 

2.  approximately  solve  the  related  problem  (4.3.2)  on  ^ 

by  calling  Outer_l_on_Mj (2);  then  make  the  correction  (4.3.3). 


The  computational  cost  of  this  multi-level  iteration,  which 
will  be  analyzed  in  the  next  section,  depends  on  the  dimensions  of  the 


spaces  {M 


In  order  to  have  a  cost  per  iteration  proportional  to 


r . . j>l . .  ^  * 

the  number  of  unknowns,  the  cost  of  the  smoothing  inner  iteration  on 
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M.  must  dominate  the  cost  of  recursively  solving  the  related  problems  on 
J0 

the  coarser  spaces.  This  will  only  be  true  if  the  dimensions  of  the 

Till00 

spaces  IM  }  ^  grow  fast  enough.  Since  these  dimensions  will  grow  quite 

1  00 

slowly  for  many  families  of  locally  refined  grids  th .  presents 

a  problem. 

To  ameliorate  this  difficulty,  we  consider  now  a  second 
approach  to  solving  the  related  problem  (4.3.2)  in  the  multi-level 
iteration  (4. 3. l)-(4. 3. 3) .  This  second  approach  is  based  on  a  slightly 
sharper  version  of  theorem  4.1.  First  we  require  a  simple  algebraic  lemma. 

Lemma  4.2.  Let  x,y  be  constants  in  (0,  1).  If  m  is  given  by 


m  = 


Lz-JL 

xy 


then 


Proof.  We  have 


(1  -  x)“  <  y  • 


. .  .m+1 

1  j-x)  /j--  =  1  +  (1  -  x)  +...+  (1  -  x)r 


>  (m  +  1)  (1  -  x)  , 


so 


Then 


so 


(1  -  x)m+1  -  1  <  -x(m  +  1)(1  -  x)m  . 


1  >  (1  -  x)^1  +  x(m  +  1 )  ( 1  -  x)m  =  (1  +  mx)(l  -  x)m  , 


(1  -  x)m  <  (1  +  mx)  <  (1  +  — m  y  *  □ 


Using  this  lemma,  we  can  show  the  following  theorem,  which  for 


our  purposes  improves  on  theorem  4.1,  All  notation  is  as  before. 


Theorem  4.3.  Let  r  £  (0,  1)  be  fixed  and  let  y  £  (0,  1)  and  j  >  2  be 
arbitrary.  Then  there  exists  a  positive  integer  m  independent  of  j 
such  that  if  the  approximate  solution  v  of  (4.3.2)  satisfies 

(4.3.8)  |||v  -  v|||  £  r  y  |||v|||  , 
then 

(4.3.9)  ||lVl  -  u'J’llI  <  Y  ll|u0  -  u0>|ll  • 

Moreover,  m  can  be  chosen  to  satisfy 

(4.3. 10)  m  £  mQ  Y 

where  m^  >  0  is  a  constant  independent  of  j  and  y. 

Proof.  The  proof  is  the  same  as  that  of  theorem  4.1,  except  for  the 
choice  of  the  parameters  m,  0^,  y  at  the  end  of  the  proof.  Instead  of 
inequality  (4.3.7)  we  have 

(4 . 3.  i  1 )  III  en)+1  III  £  (CQ0J174  +  (1  -  0.)m  +  ry)  |||e0|||  , 

2 

because  the  factor  ry  in  (4.3.8)  replaces  the  factor  y  in  (4.3.5). 
Since  we  wish  to  show 


« 


Then  we  may  choose  0^  as 


and  m  as 


0lM 

9j  ' y 


Using  lemma  4.2, 

c0e*/4  +  (1  -  6 )"  <  +  <  tt  .  r)Y , 

J  J  Q 

establishing  the  first  part  of  the  theorem.  For  the  second  part  we  have 


,  2C0  >5  -5 

ll  -  rJ  T 


completing  the  proof.  □ 


<  j  +  ,  0  -5 

—  V1  -  r>  ' 


1  u  +  (rrr)5)  Y~5  * 


Using  this  theorem,  we  can  now  analyze  a  modified  version  of  the 
multi-level  iteration  given.  This  modified  iteration  will  be  the  same 
as  (4. 3.1)— (4. 3. 3)  except  that  m,  the  number  of  smoothing  inner  iterations, 
will  now  be  made  to  depend  on  the  level,  j.  In  applying  this  modified 
iteration  to  M,  ,  j_  _>  3,  the  related  problem  in  M  will  be  solved  by 

J0  Jo_i 

doing  only  one  of  the  modified  outer  iterations  on  M  1 .  However,  to 

adequately  solve  this  related  problem  on  M  .  using  only  one  outer 

J0 

iteration,  rather  than  two  as  before,  the  number  of  smoothing  inner 
iterations  must  be  increased. 
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An  informal  statement  of  this  modified  outer  iteration 

follows.  Calling  this  procedure  with  j  =  j  causes  n  multi-level 

iterations  to  be  performed  on  M.  .  Here  m(*,  •)  is  an  integer  valued 

J0 

function  to  be  given  following  the  description  of  this  procedure. 

Procedure  0uter_2_on_M ,  (n,  j^) 

If  (j  =  1)  {solve  the  equations  on  M  directly;  then  return} 
else  {repeat  n  times{ 

1.  do  m(j ,  Jq)  smoothing  inner  iterations  on  M.  according 
to  (4.3.1). 

2.  approximately  solve  the  related  problem  (4.3.2)  on  M.  ^ 
by  calling  0uter_2_on_Mj (1,  j^),  then  make  the 
correction  (4.3.3). 

}} 

return 


It  remains  only  to  choose  the  function  m(*,  •)  required  by 
this  procedure.  Let  r  £  (0,  1)  be  arbitrary  and  let  m^  >  0  be  the 
corresponding  constant  from  inequality  (4.3.10)  of  theorem  4.3.  Also,  let 
Yq  £  (0,  1)  be  arbitrary.  Then  we  may  set 

I  -5  5(j~V 

(4.3.12)  m(j,  jQ)  =  U0  yo  r  • 

Since  the  argument  needed  to  show  that  the  multi-level  iteration  in 

0uter_2  converges  rapidly  is  somewhat  more  complex  than  that  for 

r  i 00 

Outer  1,  we  state  this  as  a  theorem.  Let  lu„  }  ,  be  successive 

—  0,n  n=l 

iterates  produced  by  the  multi-level  iteration  in  0uter_2.  Then  we  have 
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Theorem  4.4.  Let  2  be  arbitrary.  Then  with  m(#,  •)  given  by 


.  .  w 

(4.3.12)  the  iterates  iu-  )  .  satisfy 

0,n  n=l 


'  U0,n+1  "  u 


(j0) 


c  u„  -  u 

-  '0  111  0,n 


(j0^> 


for  all  n  >  1. 


Proof.  In  order  for  each  outer  iteration  on  2,  to  reduce  the 

energy  norm  of  the  convergence  error  by  the  factor  y  ,  according  to 

theorem  4.3,  it  suffices  to  do  m  of  the  inner  iterations  (4.3.1)  on 

M  for  some 
J0 


m  <  m0  Y0  , 

and  then  approximately  solve  the  related  problem  (4.3.2)  with  relative 
accuracy  rYQ.  Examination  of  (4.3.11)  shows  that  doing  extra  inner 
iterations  never  hurts  so  the  choice 


' j n >  3n>  =  l" 


m  =  V  =  lroo  Yo 

is  fine.  Thus,  for  the  case  where  =  2  and  the  related  problem  on 

M  .  is  solved  directly  the  proof  is  done.  Otherwise,  we  solve  the 
J0 

related  problem  on  M  with  one  outer  iteration.  Choosing  Y  =  rYg 

in  theorem  4.3  it  suffices  to  do  m  inner  iterations  on  M  .  for  some 

J0 


m  <  rag (tYq) 


-5 


and  then  solve  the  related  problem  on  M.  -  with  relative  accuracy 
2  j°' 


r  yq-  Doing 


m  -  m(j0~l,  jQ) 


L0  y'5  r_5| 


inner  Iterations  is  certainly  sufficient,  so  the  proof  is  done  for 
j  =  3.  Continuing  the  induction  in  the  obvious  way  completes  the  proof.  □ 


C 

'i 

M 


4 
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The  computational  cost  of  the  multi-level  iterations  described 
in  this  section  and  construction  of  an  attractive  overall  algorithm 
for  obtaining  solutions  to  the  elliptic  problem  (2.1.3)  will  be 
considered  in  the  following  section. 


4.4.  Complexity 

In  this  section  an  estimate  of  the  computational  cost  of  the 

multi-level  iterations  described  in  the  last  section  is  given,  and  these 

iterations  are  then  used  as  building  blocks  for  computationally  attractive 

algorithms  for  solving  finite  element  equations  on  locally  refined  grids. 

Under  certain  restrictions  these  algorithms  will  be  shown  to  be  of 

optimal  order,  producing  a  good  solution  to  the  finite  element  equations 

for  M^,  k  1 ,  in  (1(11^)  operations,  where  is  the  dimension  of  the 

finite  element  space  M^.  This  complexity  analysis  is  somewhat  harder 

than  the  corresponding  analysis  in  chapter  three,  since  the  complexity 

of  the  multi-level  iterations  in  the  last  section  depends  on  the 

dimensions  of  the  finite  element  spaces  {M . ,  used.  This  complication 

J  J  =  1 

was  absent  from  the  last  chapter,  since  for  the  family  of 

quasi-uniform  spaces  there,  the  dimensions  of  the  spaces  were  essentially 

determined  by  the  dimension  of  the  coarsest  space  and  the  mesh  ratio  p . 

A  second  minor  difficulty  also  arises  in  proving  the  optimal 

r  '•-!  divergence  of  the  algorithms  here.  A  multi-level  algorithm  is 

. 'red  to  lie  of  optimal  order  if,  when  applied  to  any  finite  element 

~(k) 

i-  •'!  ,,  it  produces  a  good  solution  uv  in  (1(N,  )  operations. 

11  =  1  k 

.  "t  good  solution"  means,  an  error  bound  applicable  to 

■  t  ills  is  needed.  We  now  proceed  to  derive  such  a  bound. 


r  hat  t  hi"  error  bound  derived  not  be  overly 


pessimistic.  If  it  were,  the  results  here  would  be  correspondingly 
weakened.  To  have  shown  a  multi-level  method  produces  a  solution 
containing  a  convergence  error  smaller  than  the  expected  truncation 
error  means  very  little  if  one  expects  far  more  truncation  error  than 
actually  results. 

Only  a  simple  error  estimate  is  required  here.  Such  an 

estimate  can  be  derived  directly  from  the  Bramble-Hilbert  lemma.  First 

we  need  the  notion  of  a  Sobolev  seminorm.  For  £  a  non-negative  integer, 

the  seminorm  |  •  |  is  the  same  as  the  Sobolev  norm  ||  •  ||  .  except  that  it 
H  IT 

contains  only  derivatives  of  order  £.,  not  those  of  lower  order.  That 

is ,  for  u  S  , 

|u|  *  -  <  i  l»A.H22>1/2  . 

H  (fi)  j  ot  |  =£  L 

Just  a  for  the  Sobolev  norm,  these  seminorms  can  be  defined  for 
non-integral  I  in  a  more  complicated  way. 

Now  let  M  be  a  finite  element  space  of  Lagrangian  elements, 
based  on  a  triangulation  T  of  ft  and  satisfying  conditions  (4.2.1)-(4. 2.5) . 
Suppose  the  element  trial  functions  are  polynomials  of  degree  k  -  1. 

Then  M  is  said  to  be  a  finite  element  space  of  degree  k.  For 

l 

u  £  H  (ft) ,  1  <  £  <  k,  the  nodal  interpolant  u  £  M  of  u  is  well  defined 

h  *—  I 

and  we  have 


Lemma  4.5.  For  any  triangle  T  £  T  and  for  s  =  0,  1 

lu  -  uxi  s  £  °s  4”s  M  a 

1  H  (T)  1  H  (T) 


where  d^  is  the  diameter  of  T. 
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The  variable  l  here  need  not  be  an  integer.  This  lemma  is  a  special 
case  of  the  Bramble-Hilbert  lemma,  which  may  be  found  in,  for  example, 
Strong  and  Fix  or  Oden  and  Reddy. 

.  >00 

Applying  this  lemma  to  the  spaces  in  a  family  {M  }  ^  °f 

locally  refined  finite  element  spaces  satisfying  the  assumptions  of 

section  4.2  yields  the  required  error  bound.  Assume  that  the 

differential  equation  (2.1.3)  is  H  regular  for  ctG  (0,  1),  that 

^  |  ^ 

is,  that  its  solutions  lie  in  H  .  Then  the  following  error  estimate 
holds  on  M . ,  i  >  1 . 

j  - 


Theorem  4.6.  Let  uu  be  the  finite  element  solution  of  (2.1.3)  in 


M. .  Then 
1 

(4.4.1) 


|u  -  u(j)m  <  c  ( i  d,^a  lul21+a  y 

TGT  1  H.  (T) 


where  c^  >  0  is  independent  of  u,  j,  and  the  family  of  locally  refined 


r  1  t  00 

spaces  {M . } .  - . 

ii=l 


Proof .  Let  be  the  interpolant  of  u  in  .  Then 

III  U  -  u^HI  =  min  |j|u  -  v||j  <.  |||u  -  u^ 
v^M 


Since  the  energy  and  H  norms  are  equivalent. 


(4.4.2) 


N  -  «lJ)lll  i  •  IN  -  U!J)II  !  • 

H 


Applying  the  above  lemma  on  each  triangle  in  , 


(j)i|2  <  2  ..  .2(140)  |  1 2 

'[  I*  2  -  0  1  d7  |u!  140 

1  1.(11)  UTei.  H  (T) 
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Thus , 


4J)i2i  i  1  d2ct 

1  H  (£1)  1  TSt.  T 

J 


I  (j)|| 2  2  2a 

1“  ~  u II  1  i  c  ^  dT 
h  (fl)  ieT.  1 

j 


2  2 , 


H1+Ct(T) 


H1+0t(T) 


for  c  =  max{a0,  a^},  assuming  dT  £  1  for  all  triangles  T  €  x  , 
combining  this  with  (4.4.2)  the  result  follows.  □ 


Now 


Since  this  bound  contains  d^,  in  effect  the  local  mesh  parameter, 

it  is  much  tighter  for  locally  refined  grids  than  bounds  in  terms  of  the 

maximum  element  diameter.  Moreover,  the  exponent  of  d^  here  is  optimal. 

For  sufficiently  smooth  u  the  reverse  inequality,  bounding  the  error 

from  below,  may  be  expected  to  hold  although  this  is  difficult  to  show. 

Establishing  that  the  multi-level  methods  of  this  section  produce 

solutions  whose  convergence  error  is  comparable  to  the  truncation  error 

expected  from  this  error  bound  may  be  regarded  as  more  than  adequate  for 

most  purposes.  In  particular,  for  problems  with  singular  solutions 

the  improved  convergence  of  finite  element  spaces  using  the  proper 

locally  refined  grids  may  be  readily  shown  using  this  bound.  The 

multi-level  algorithms  here  will  inherit  this  improved  convergence  to 

singular  solutions  on  properly  refined  grids. 

With  this  error  bound  established  we  are  ready  to  turn  to  the 

design  of  multi-level  algorithms  for  solving  the  elliptic  problem  (2.1.3) 

on  locally  refined  grids.  Let  be  a  family  of  locally  refined 

finite  element  spaces  satisfying  the  assumptions  of  section  4.2.  The 

~  00 

algorithms  we  wish  to  consider  will  produce  a  good  solution  u  in  any 
number  of  this  family,  k  1.  The  first  step  in  these  algorithms  is 


I 


4 


81 


to  solve  the  finite  element  equations  for  directly  to  obtain  an 
approximate  solution  to  the  problem  (2,1.3).  After  this,  solutions 

£  M  ,  2  £  j  <_  k,  will  be  produced,  each  a  good  approximate  solution 
of  (2.1.3).  Each  solution  2  <  j  <  k,  is  obtained  by  beginning 

with  the  previous  solution  u^  ^  £  A).  ^ ,  considering  it  instead  as  a 
function  in  M  ,  and  then  carrying  out  a  sequence  of  multi-level  iterations 
to  obtain  the  improved  solution  u^  €  M . .  This  multi-level  iteration 


J 

may  be  taken  as  either  of  the  algorithms  of  the  last  section. 

Two  questions  arise  immediately:  what  conditions  must  the 

~  (k) 

above  process  satisfy  in  order  to  obtain  an  accurate  solution  u  , 

and  what  determines  the  computational  complexity  of  this  process?  With 

the  error  bound  (4.4.1)  already  proven,  the  first  question  is  easier. 

~  (k) 

We  wish  to  obtain  a  solution  u  which  contains  only  convergence  error 
comparable  to  the  expected  truncation  error.  That  is,  we  want  the 
convergence  error  to  satisfy 


lllu00  -  „<k)|||  <  c  <  E  df  |u| 

k 

In  order  to  have  an  algorithm  whose  solution 
inequality,  it  is  convenient  to  ask  that  the 
u(j),  l  <_  j  £  k  -  1,  also  satisfy  this  type  o 
for  1  <  j  <_  k, 

(4.4.3)  |||a(j)  -  u(j)|||  <  c  (  Z  d2a  |  u  j 

TGx.  1 
] 


,1/2 


H1+a(T) 


~(k)  .. 

u  satisfies  this 


intermediate  solutions 
f  inequaltiy.  That  is. 


,1/2 


H1+a(T)> 


A  simple  condition  insuring  this  is  given  in  the  following  theorem. 
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dT  <  2dT,  ,  T'  G  ST 

I u  1 2H-ot  =  Z  lu|21+a 

H  (T)  T’GS  HI+a(T')  . 

Applying  these  to  (4.4.5)  yields 

Hl=a>  -U<j+1>llllc0(  E  df  |u)2  )1'2 

T^Xj+1  H  (T) 

+  2c  (  Z  Z  (2d  )2a  |u|2  )1/2 

U  TGT.  T'GSt  h  (t') 

=  c  (1  +  2 (2)a)  (  Z  d2a  J u I  )1/2  . 

°  T6Tj+1  T  H1+a(T’) 

Finally,  making  use  of  the  hypothesis  (4.4.4)  the  induction  is  complete.  □ 


The  dependence  of  this  result  on  the  regularity  of  the  problem 
will  affect  the  complexity  of  the  multi-level  algorithm  here  in  an 
interesting  way.  Consider  the  construction  of  such  algorithms  based  on 
this  theorem.  As  described  above,  after  obtaining  u^^  by  solving  the 
equations  for  directly,  approximate  solutions  u^  £  M.,  2  <  j  <  k, 
are  obtained  by  iteratively  improving  the  previous  solution  u  ^  ^  6  M 
The  amount  of  iterative  improvement  required  is  given  by  (4.4.5).  If  the 
problem  is  H  regular,  we  must  reduce  the  convergence  error  by  the 
factor  (1  +  2^  *)  ^  in  going  from  ^  to  u^\  The  interesting  point 

here  is  that  more  error  reduction  is  required  for  smoother  problems. 

Since  the  regularity  did  not  enter  the  convergence  theorems  of  the  last 
section,  obtaining  a  given  amount  of  reduction  of  the  convergence  error 
is  apparently  no  easier  for  problems  with  smooth  solutions  than  for 
problems  with  weak  regularity.  Thus,  the  computational  cost  per  element 


of  the  algorithms  here  will  be  greater  for  problems  with  smooth  solutions. 
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Of  course,  there  is  no  real  dilemma  here.  The  algorithms 

here  produce  a  solution  satisfying  (4.4.5).  Since  the  finite  element 

(k) 

solution  u  satisfies  (4.4.1),  the  approximate  solution  produced  by 
the  multi-level  algorithm  must  satisfy 


|u(k)  -  u|||  <2c  (  E  d2a  |u|2 

°T€xk  T  H14<X(T) 


,1/2  , 


a  much  stronger  inequality  for  large  a.  Thus,  even  though  one  must  work 
harder  per  element  for  smooth  problems,  it  will  still  be  cheaper  to 
produce  a  good  answer. 

Now  we  are  ready  to  consider  the  design  of  multi-level 
algorithms  based  on  theorem  4.  7  and  the  multi-level  iterations  of 
section  4.3.  The  multi-level  iteration  used  may  be  that  described  in 
either  procedure  Outer_l  or  0uter_2.  In  either  case  the  dimensions  of 
the  finite  element  space  1M.K will  be  of  central  importance.  If  N^ 
is  the  dimension  of  M  ,  j  ^  1,  optimal  complexity  bounds  can  be  shown 
under  the  assumption  that  the  dimensions  {N ^ 1 j _ ^  grow  at  least 
geometrically.  That  is, 


(4.4.6) 


Vi  i  ,u>. 


for  fixed  6  >  1.  Actually,  this  can  be  weakened  slightly.  It  is  enough 

r  l°° 

that  we  have  a  sequence  of  numbers  growing  at  least  geometrically, 

XJ+1  ±  6xj  ’  j  ^  1  * 

for  (5  >  1,  and  that  and  x^  are  comparable, 

—  x  <  N .  <  c  x 
c  j  ~  J  -  j 

for  fixed  c  >  0.  For  simplicity,  assume  that  (4.4.6)  holds.  Using  either 


of  the  multi-level  iterations  of  section  4.3,  the  storage  required  to 
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approximately  solve  the  finite  element  equations  for  k  >  1,  can  be 
estimated  as  follows.  Let  c^  be  the  storage  required  for  direct 
solutions  of  the  equations  for  .  Then  the  total  storage  required 
is  proportional  to 


k.  •  t 

Z  N.  +  c  <  N,  Z  +  c, 

j=2  J  k  j=2  1 


<  Nk 


1  “  6 


1  +  C1 


Since  c^  is  fixed,  for  any  6  >  1  the  storage  required  is  proportional 

to  N.  • 
k 

Next  consider  the  use  of  the  multi-level  iteration  described 
by  procedure  Outer_l.  According  to  theorem  4.1,  when  this  multi-level 
iteration  is  applied  to  ,  j  >  1,  the  convergence  error  is  reduced  by 
a  factor  Y  at  every  iteration,  where  y  <  1  is  independent  of  j. 

Actuall  the  proof  of  theorem  4.1  shows  Y  can  be  chosen  arbitrarily, 
though  this  is  unimportant,  since  for  any  y  <  1  we  can  find  v  _>  1 
such  that 


v  ^  ,  -1-fa .-1 

Y  1  (1  +  2  ) 

Let  y  <  1  and  an  integer  v  satisfying  this  inequality  be  fixed.  Then 

~(k) 

according  to  theorem  4.7,  a  solution  u  in  M^,  satisfying  (4.4.5)  can 
be  produced  by  carrying  out  v  of  the  multi-level  iterations  described 
by  Outer_l  on  each  space  ,  2  <  j  £  k.  For  each  ,  2  j  k,  the 
approximate  solution  produced  in  the  next  coarser  space  is  used 

as  a  starting  value  for  this  iteration.  To  estimate  the  complexity  of 
this  procedure,  note  that  in  applying  Outer_l  to  M  ,  j  _>  2,  O(N^) 
operations  are  required  plus  the  cost  of  recursively  solving  the  related 


« 


» 
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problem  in  j.  This  related  problem  is  solved  in  Outer_l  applied  to 
M.,  by  recursively  applying  Outer  1  twice  to  the  equations  for  ^ ,  if 
j  -  1  >  1.  Continuing  the  recursion  one  sees  that  the  total  cost  in 
doing  V  multi-level  iterations  of  this  tvpe  on  the  equations  for  is 
proportional  to 

v(  £  2j_i  N.  +  2J_1  c.)  <  VN.(  Z  (4)j-1  +  (|)j_1  c„)  , 


1=2 


”2  - 


i=2 


where  c9  is  the  cost  of  solving  the  linear  system  for  Since  is 

a  fixed  constant,  the  computation  time  T(M^)  t0  do  V  of  these  multi¬ 
level  iterations  on  M.  must  satisfy 


T(M. )  < 


<,(V  ■ 

0(N  log  (N  )) 

log  2 

108  S  , 


6  >  2 
6  =  2 

1  <  6  <  2  . 


Then  following  the  procedure  described  above  where  the  approximate  solution 

~  00 

u  in  Mk  is  obtained  by  successively  obtaining  good  solutions  in 

~(k) 

M  ,  1  <  j  £  k,  the  total  computation  time  T  to  obtain  u  must  satisfy 


<V  • 


6  >  2 


T  =  Z  T(M. )  <  J  (N.  log  (N  ))  ,6=2 
J-l  J  '  k 

log  2 


,(»k  108  s>  • 


1  <  6  <  2  . 

For  many  problems  a  bound  like  this  is  quite  adequate.  For 
example,  in  doing  local  refinement  along  a  plane  singularity  or  boundery 

.  OO 

in  a  three  dimensional  problem,  the  dimensions  would  at  least 

quadruple  at  each  level.  Since  the  two  dimensional  results  here 
readily  extend  to  three  dimensional  problems,  this  example  is  important. 
In  other  problems,  for  example,  refining  along  a  line  sigularity  in  two 


' 
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or  three  dimensions,  the  dimensions  typically  only  double  at  each  level, 
yielding  the  N  log  (N)  bound,  Stil1  other  problems  arise  for  which  the 
growth  of  these  dimensions  is  even  slower.  In  refining  around  point 
singularities,  these  dimensions  will  usually  not  grow  geometrically, 
and  the  results  of  this  section  are  of  little  use.  This  does  not 
necessarily  mean  that  multi-level  methods  should  be  avoided  in  these 
cases,  only  that  theoretical  results  showing  their  effectiveness  have 
not  yet  been  given. 

When  the  dimensions  of  the  spaces  M.  grow  geometrically 
according  to  (4. A. 7),  but  the  rate  of  growth,  6,  is  not  greater  than 
2,  the  complexity  results  above  can  be  improved  by  replacing  the 
multi-level  iteration  given  by  Outer_l  with  that  given  by  0uter_2. 

In  equation  (4.3.12),  governing  the  number  of  inner  iterations  used  on 
each  level  in  0uter_2,  we  are  free  to  choose  both  y^  and  r.  For 
convenience,  y^  may  be  chosen  as 

y0  -  (i  +  21+ct)-1  . 


With  this  choice  in  applying  0uter_2  to  the  equation  for  M.  , 

^  0(  j  ) 

2  <_  Jq  <_  k,  only  one  iteration  of  0uter_2  is  needed  to  get  u  0 
satisfying  (4.4.4)  from  u  ^  .  The  other  free  parameter,  r,  in 

(4.3.12)  will  be  chosen  as 


=  6 


Then  for  any  6  >  1,  re  (0,  1),  so  theorem  4.4  shows  that  applying  the 


multi-level  iteration  described  by  0uter_2  to  M  ,  2  _<  j  <_  k,  reduces 

J0  0 


the  convergence  error  by  y^  per  iteration  as  required.  From  equation 
(4.3.12),  governing  the  number  of  inner  iterations  on  each  level,  we  can 
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Ji0)  JV1) 

estimate  the  computation  time  to  produce  u  from  u  ,  ^  2. 

This  time,  T(M.  )  is  proportional  to 
J0 

J0  -5(j  -j)  J0  -5  Vj 

E  r  N.  +  c0<  N.  -E  (—  )  +  c 

J-2  2  2  ”  J0  j-2  «  2 


VJ  A  *  —  ’ 

<  N.  E  (6  b)J0  J  +  c 
J0  j=2  ] 


<  N.  (1-6  V1  +  c  . 
“  J0  2 

<V 


Thus,  the  amount  of  computation  required  to  produce  u  from  u 


(V1) 


j  >  2,  is  t? (N .  ).  Then,  since 
‘  J0 


K  -1  -1 
E  N.  <  N  (1  -  6  ) 

J„-2  J0 

for  any  6  >  1,  the  total  computation  time  is  O(N^). 

This  is  an  optimal  bound,  showing  0(N)  complexity  under  the 

.  OO 

weaker  assumption  that  the  dimensions  iN  grow  geometrically,  with 

j  3  ^ 

any  growth  rate  6  >  1.  Unfortunately,  the  dependence  of  this  ()(N)  bound 
on  r;  is  quite  severe,  though  this  could  be  improved  by  choosing  the 
various  constants  differently,  and  by  using  a  more  effective  inner 
iteration.  The  hidden  constant  in  the  0(N)  bound  just  given  contains 
not  only  the  factors  (1-6  ^  and  (1  -  6  ^  seen  above,  but 


also  the  constant  m.  of  (A. 3. 12),  which  depends  indirectly  on  6.  By 
U  ) 

the  proof  of  theorem  A. 3  it  can  be  seen  that  m^  goes  as 


89 


« 


_  1 

(1  -  r)"5  =  (1  -  5  V5 

Thus,  including  the  dependence  on  6,  the  complexity  bound  is 

_  1 

0((1  -  6  6)'6  (1  -  6"1)'1  N)  . 

While  the  factor  (1-6  (1  -  6  ^  already  exceeds  10^  for 

6=2,  in  practice  one  would  expect  the  dependence  on  6  to  be  far 
less  severe  than  this,  although  the  appropriate  numerical  experiments 
have  not  yet  been  carried  out. 

This  completes  the  analysis  of  the  algebraic  aspects  of  the 
multi-level  algorithms  of  this  chapter.  It  remains  to  prove  the 
approximation  result  (4.3.4),  on  which  the  complexity  results  just 
given  are  based. 


4.5.  Interpolation 

This  section  begins  the  analysis  of  the  approximation  properties 
of  locally  refined  finite  element  spaces  needed  to  justify  the  convergence 
and  complexity  results  of  the  last  two  sections.  For  convenience  this 
analysis  is  split  into  two  parts.  The  first  of  these,  given  in  this 

section,  deals  with  properties  of  the  interpolation  mapping  from  Mj  to 

00 

j  >  2,  for  Mj ,  Mj  ^  members  of  the  family  {M  }  ^  of  locally 
refined  spaces.  The  second  part  of  this  analysis,  given  in  the 
following  section  will  use  the  results  of  this  section  to  obtain  the 
required  approximation  theorem. 

While  the  analy  is  in  the  next  section  is  completely  general, 
the  analysis  here  is  specific  to  the  triangular  elements  and  grids 
based  on al  igned  triangulations  being  considered.  However,  there  is  no 


reason  to  believe  that  similar  results  could  not  be  shown  for  more  general 
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grids  and  for  most  common  elements  Including  triangular  elements  with 
higher  continuity,  square  or  quadrilateral  elements,  and  elements  with 
curved  edges.  The  results  here  were  not  proven  in  greater  generality 
primarily  because  it  did  not  seem  to  aid  understanding  of  the  multi-level 
algorithms  here  and  because  the  notations  and  arguments  required  would 
have  been  burdensome  to  both  author  and  reader.  Proof  of  results 
analogous  to  those  in  this  section  for  most  reasonable  finite  element 
spaces  is  not  difficult. 

The  work  in  this  section  is  closely  related  to  the  work  Bank 
and  Dupont  did  in  analyzing  their  two-level  scheme  discussed  briefly  in 
section  4.1.  Much  of  the  notation  here  and  one  of  the  lemmas,  lemma  4.10, 
are  taken  directly  from  their  paper.  The  difference  here  is  mainly  in 
emphasis.  While  lemma  4.10  plays  a  fairly  small  role  here,  It  is 
central  to  the  convergence  proof  for  the  two-level  scheme  in  their 

paper.  Lemma  4.11,  which  is  not  from  their  paper,  is  really  much  more 

central  here. 

The  reason  for  considering  first  the  interpolation  mapping  is 
the  usual  one:  while  projection  with  respect  to  the  energy  inner  product 
is  global,  interpolation  gives  a  local  projection  more  easily  analyzed. 

*  t.  t  he  two  mappings  are  sufficiently  al  ike  that  properties  of  one  can  be 
derived  from  the  other.  This  is  a  common  point  of  view,  hut  the  way  in 

which  properties  of  the  elliptic  projection  are  derived  here  from  the 

interpolation  projection  is  unusual.  Instead  of  bounding  the  approxi¬ 
mation  error  in  the  interpolate  and  using  this  as  a  bound  on  the 
approximation  error  in  the  elliptic  projection,  we  hi  gin  hv  establishing 
properties  of  the  null  space  of  the  interpolation  mapping.  From  these 


« 


>D-AO08  065  ILLINOIS  UN IV  AT  UR BAN A -CHAMPAIGN  DEPT  OF  COMPUTER  SCIENCE  F/g  lt/S 

RAPID  SOLUTION  OF  FINITE  ELEMENT  EOUATIONS  ON  LOCALLY  REFINED  g— ETCIU) 
MAY  GO  JR  VAN  ROSENDALE  AFOSR-75-2G54 

UNCLASSIFIED  UlUCDCS-R-GO-lOai  AFOSR-TR-80-0541  ML 


91 


properties,  in  the  next  section  we  will  derive  properties  of  the  null 
space  of  the  energy  projection.  Finally  these  lead  to  the  required 
approximation  result. 

This  proof  may  seem  circuitous,  but  the  author  was  unable  to 
establish  an  approximation  result  such  as  (4.3.4)  in  any  other  way. 

The  problem  encountered  with  the  usual  finite  element  approach  is  not 
hard  to  see.  Since  elliptic  projection  is  global,  an  error  estimate, 
such  as  theorem  4.6,  for  a  locally  refined  space  M,  must  involve  the 


maximum  element  diameter  h  .  However,  the  eigenvalue  bounds  involve 

max 

h  .  .  Thus,  one  will  be  able  to  prove  inequalities  such  as  (4.3.4)  but 
mm  h 

ITIclX 

the  constant  c  occuring  will  contain  some  power  of  the  factor  - -  . 

min 

This  factor  becomes  unbounded  for  the  locally  refined  grids  being 
considered,  so  such  a  result  would  be  too  weak  to  prove  (?(N) 
complexity  bounds  like  those  in  the  last  section. 

Since  the  finite  element  spaces  considered  in  this  chapter 
are  based  on  C®  Lagrangian  elements,  the  interpolation  mapping  from 

DO 

to  M.  ,,  j  >  2,  with  M.,  M.  ,  members  of  {M.}.  ,,  is  automatically  well 
J-l  ~  j  j-1  J  J=1 

defined.  For  any  u  £  ,  u  is  evaluated  at  the  points  corresponding  to 

the  nodal  parameters  of  Then  there  is  a  unique  u^  £  attaining 

the  same  values  at  these  points.  One  of  the  problems  encountered  when 
using  elements  with  higher  continuity  is  that  this  interpolation  may  not 
be  well  defined.  There  are,  for  example,  elements  using  first 
derivatives  as  nodal  parameters  and  almost  all  elements  use  second 
derivatives.  In  these  cases  care  must  be  taken  that  the  nodes  of 
lie  at  points  where  the  corresponding  values  of  functions  in  M  are 
well  defined.  In  those  rare  cases  where  nodal  interpolation  is  not  well 


» 


defined,  more  complex  interpolation  projections  based  on  various  local 
averaging  techniques  must  be  used. 


Now  let  M_.  ,  j  2 ,  denote  the  null  space  of  the  interpolation 

A 

from  M.  to  M.  M.  is  the  subspace  of  M.  consisting  of  functions 
j  j-1  J  ■  j 

vanishing  at  the  nodes  of  Much  of  this  section  and  part  of  the 

next  are  devoted  to  establishing  properties  of  these  null  spaces.  For 
now  we  note  that  $  yields  a  direct  sum  decomposition  of  M : 


■1 .  =  M .  ©  M  • 

J  J  w  j-1 


This  is  easy  to  show,  for  example,  by  dimensionality  arguments. 

Besides  the  direct  method  of  defining  the  interpolate 
Uj  €  °f  a  function  u  G  ,  there  is  also  a  more  abstract  way  which 

we  now  consider. 

The  interpolation  mapping  from  to  considered  here 

is  completely  local  in  the  sense  that  the  value  of  the  interpolate 
“i s  Vi  on  a  triangle  T  G  Tj-i  depends  only  on  the  value  of  u  6 
on  T.  For  this  reason  the  interpolation  mapping  from  to 
decomposes  into  separate  mappings  from  M  |^,  to  ^  for  each 

T  £  Tj  ^ .  Rather  than  considering  separately  each  of  the  function 
spaces  I y  and  M  as  T  ranges  over  M  ^  we  can  consider  only  a 

few  canonical  function  spaces  on  T^  and  the  interpolation  mappings 
between  these  canonical  spaces.  These  canonical  spaces  on  T^  can  then 
be  transformed  into  the  spaces  M ^  | and  through  the  affine 

transformation  F^,. 

As  in  section  4.2  let  L  be  the  element  trial  space  on  the 

reference  triangle  T„.  We  must  consider  three  different  refinements 

of  Td:  regular  refinement,  green  refinement,  and  the  null  refinement, 

K 
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in  which  TR  is  not  altered.  There  are  three  possible  green  refinements 
but  we  need  only  one,  say  the  one  with  the  new  edge  leading  to  the  top 
vertex. 


Regular  Green  Null 

Figure  16.  Refinements  of  T 

K 


Since  each  of  these  refinements  of  T  is  a  triangulation, 

K 

finite  element  spaces  on  each  of  them  can  be  generated  in  the  usual  way 

from  the  reference  element  trial  space  L. 

Let  K  be  the  finite  element  space  on  T„  generated  in  this 
K  K 

way  from  the  regular  refinement  of  T  and  let  K  be  the  finite  element 

K  Cj 

space  generated  by  this  green  refinement.  Null  refinement  just 

reproduces  L.  All  three  of  these  spaces  generated  by  refinement,  K  , 

R 

Kq,  and  L  contain  L  as  a  subspace.  Moreover,  there  are  interpolation 
mappings  IR,  IG,  and  1^  from  KR,  K G,  and  L,  respectively  to  L.  The 
last  is  of  course  just  the  identity  on  L.  Just  as  for  the  interpolation 


mapping  from  M.  to  M,  , ,  we  denote  the  kernels  of  ID  and  I,,  by  K  and 
J  j-l  K  b  K 

Kf, ,  respectively.  The  kernel  of  1^  is  the  trivial  function  space  (0). 

Now  let  T  E  T.  ^  be  arbitrary.  In  going  from  to  Tj  the 

triangle  T  is  regular,  green,  or  null  refined.  Consider  the  same  type 


of  refinement  of  T  .  Then  there  is  an  affine  transformation 
R 

F  :  T  -*■  T  carrying  the  subtriangles  of  T„  generated  by  this  refinement 
IK  K 

to  the  subtriangles  created  in  refining  T.  When  T  is  regular  or  null 
refined,  any  of  the  six  affine  transformations  have  this  property.  When 


T  is  green  refined  two  of  the  six  affine  transformations  carry  the  green 
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subtriangles  of  TR  to  those  of  T.  The  other  four  affine  transformations 

correspond  to  the  other  possible  green  refinements  of  T. 

To  see  the  importance  of  all  this,  let  IT  be  the  interpolation 

mapping  from  | ^  to  M  T  and  let  F^,  he  the  dual  of  FT,  carrying 

functions  on  T  to  functions  on  T  , 

R 

F^,:  H1  (T)  H1(Tr) 


by 


Fj(u)  =  u  o  FT  . 

For  each  of  the  three  types  of  refinement  of  T  we  get  a  commutative 
diagram.  For  regular  refinement  we  have 


t1* 


F’ 

T 


V.'t 

t1! 


For  green  refinement  we  have 


F' 

T 


M 


j  T 


r< 


F' 

T 


j  - 1  T 


j  T 


And  finally  for  null  refinement  we  have 
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F' 

T 


th 


Vi't 

tXT 


F' 

T 


"j't 


That  these  diagrams  are  in  fact  commutative  is  a  consequence 
of  the  fact  that  F^  carries  not  only  TR  and  its  subtriangles  into  T  and 
its  subtriangles ,  but  carries  nodes  of  TR  and  its  subtriangles  into 
corresponding  nodes  of  T  and  its  suhtriangles.  We  omit  the  straight¬ 
forward  but  messy  verification. 

The  lemmas  which  follow  just  elucidate  the  properties  of 
these  three  commutative  diagrams  and  the  various  function  spaces  involved. 
The  first  result  here  shows  that  functions  in  L  cannot  well  approximate 

nonzero  functions  in  K  or  K  ,  the  kernels  of  I  and  I  .  The  kernel, 

K  U  K  G 

{o},  of  1^  contains  no  nonzero  functions  so  the  question  does  not  arise 
there.  The  functions  in  KR  and  vanish  at  the  nodes  of  TR  while 
nonzero  functions  in  L  cannot  vanish  at  all  nodes  of  TR.  The  results 
follow  from  this  and  the  finite  dimensionality  of  all  spaces  involved. 

For  convenience  define  the  interpolation 

vL:  C°(Tr)  +  L  . 


The  mappings  IR,  ,  and  1^  described  above  are  just  restrictions  of 
to  the  subspaces  KR,  ,  and  L  of  C^(TR),  respectively. 


A  A 

Lemma  A. 8.  If  w  £  K„  or  w  6  ft  ,  for  any  v  €  L 

_  ■  K 
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(4.5.1a) 


(4.5.1b) 


|w  -  v||  2  >  c  ||v|| 

L  (Tr)  l/(TR) 


w  -  v|  L  >  c  |w|  1 


H  <V 


H  <V 


for  c  >  0  Independent  of  v  and  w. 

Proof.  Let  K  =  +  K„.  Then  K  is  contained  in  the  null  space  of  rr , . 

-  K  L 

Since  u.  is  the  identity  mapping  on  L, 

(4.5.2)  K  n  1  =  ■$  . 

Now  consider  the  functional 

f(w,  v)  =  || w  -  v||  , 

L  (Tr) 

defined  on  K  x  L.  Let  C  be  the  compact  subset 

{(w,  v):  1 1 w ||  =  |j v ||  =1} 

L  (TR)  L  (Tr) 

A 

of  K  x  L.  Since  f  is  continuous  it  attains  its  intimum  of  some  (wQ,  vQ)  6  C, 
but  by  (4.5.2)  we  have 


'  V°l2<v i  c  >  0 ' 


or 


lw0  '  V0H  2  -  c  2 

°  l  (Tr)  l/(TR) 


Since  f(w,  v)  is  minimal  at  (w^,  v^) ,  (4.5.1a)  holds  on  all  of  C,  and 
then  on  k  x  L  by  linearity. 


For  (4.5.1b)  note  that  since  the  functions  in  K  are  continuous 

and  piecewise  C  ,  the  only  functions  u  £  K  for  which  |u|  .  =  0  are  the 

H 

constant  functions.  The  only  constant  function  in  k  is  the  function 
identically  zero,  since  k  is  contained  in  the  null  space  of  and  the 
function  values  at  the  vertices  of  TR  are  nodal  parameters  for  L.  Thus, 


« 


» 
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the  seminorm  [•|  1  is  norm  on  K.  Because  of  this,  the  verification  of 
(4.5.1b)  goes  through  exactly  the  same  as  for  (4.5.1a).  □ 


For  any  triangle  T  c  IR  ,  each  affine  transformation 
F^:  TR  -*■  T  induces  a  dual  mapping  F^,  carrying  functions  on  T  back  to 


functions  on  T. 


R- 


by 


F^,:  H1  (T)  +  H1(Tr) 


F^,(u)  =  u  o  Ft 


One  can  then  ask  for  the  norm  of  this  dual  mapping  F,J,  either  as  a 

mapping  from  H1(T)  to  H1 (TR)  or  from  L2(T)  to  L2(TR).  For  the  H1 

seminorm  the  answer  is  that  |u|  ,  and  | u  °  F  |  are  essentially 

H  (T)  T1 (Tr) 

the  same  since  the  change  in  area  over  which  | • |  1  is  computed  and 

H 

change  in  the  size  of  the  derivatives  are  opposite  effects  and  cancel 
perfectly  (in  two  dimensions).  Only  triangles  I  6  Tj-i  are  considered 
in  the  following  lemma  both  because  they  are  the  only  ones  of  interest 
here  and  because  the  requirements  of  section  4.2  force  all  triangles 
in  Tj  i  to  have  diameter  comparable  to  h^  simplifying  the  statement 
of  the  lemma. 


Lemma  4.9.  Let  T  6  j  ^  2  and  let  F  be  an  affine  transformation 

from  T  onto  T.  If  u  €  H^(T),  then  u*  =  u  0  F  £  H^(T„)  and 
K  K 

(4.5.3a)  7  M  ,  1  lu* I  ,  1  c  |u|  . 

H1 (T)  H  (TR)  H  (T> 


t 


9 


(4.5.3b)  ^  ||u||  <  h  ||u*||  , 

I*  (T)  J  L  (Tr) 


<  c  u 


1  2 

L  (T) 


for  fixed  c  >  1  independent  of  u,  j,  and  T. 

2 

Proof .  F  =  Ax  +  b  for  some  2x2  matrix  A  and  b  £  R  .  Then  as  shown 
in  Oden  and  Reddy  (1976,  eqs.  (6.167)  and  (6.168)),  for  any  s  >  0 


(4.5.4a) 


»S<V 


<  IWf  Idet  Af1/2  |u| 


HS(T) 


(4.5.4b) 


HS(T) 


<  || A’1  If  |det  A|1/2  |u*| 


hS(Tr) 


where  ] | A | [  is  the  spectral  norm.  Now  let  d  and  d  be  the  diameters  of 

1  TR 

T  and  T  ,  respectively  and  let  p„,  and  p  be  the  diameters  of  the 
inscribed  circles  of  T  and  TR.  Then  as  given  in  Oden  and  Reddy, 
lemma  6.2,  the  following  simple  bounds  on  the  norms  of  A  and  A  1  apply 


We  also  have 


lA-'lli- 


det  A  = 


area  (T) 
area  (TR) 


The  triangle  T  satisfies  the  angle  condition  (4.2.1)  uniformly  for 
T  £  Tj-l  ant*  J  1.  2.  Using  this  condition  it  can  be  shown  by  simple 
geometric  arguments  that 


—  c  pm 

c  T  -  mT 


—  d  <  p_, 
c  T_  —  mT 
R  R 

^  <  area(T)  <  c  d2 

7  d2  1  area(TR)  <  c  d2  , 

R  R 

for  c  >  1  independent  of  T  £  t,  n  and  j  ^  2.  Then  since  d  =  1  by 

3-i  Tr 

assumption  it  follows  that 

II A H  <  c  d„ 


II  A_1 1|  £  c  dxL 
| det  A|1/2  <  c  dT 
I det  A| -1/2  <  c  dj1  . 

Combining  these  relations  with  (4.5.4)  yields  (4.5.3a)  and 


7  IMI  2  <  d  ||u*|| 

L^(T)  L  (TR) 


<  c  u 


1  2 

L  (T) 


Since 


—  h .  ,  <  d„  <  c  h .  , 
c  j-1  -  T  J-1 

by  (4.2.  ),  and  p  =  2,  inequality  (4.5.3b)  follows  and  the  proof  is 
complete.  □ 


Lemma  4.9  above  shows  how  the  norm  of  a  function  u  on  a  triangle 


T  is  affected  by  transforming  it  to  a  function  u*  ■  u  0  on  TR. 

Lemma  4.8  showed  that  functions  in  K  or  K  were  not  well  approximated 

R  v  • 

by  those  in  L.  Combining  these  two  lemmas  we  should  be  able  to  show 
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that  for  any  triangle  T,  functions  in  Z  0  F  ^  or  K  o  F  ^  are  not 

K  T  G  T 

well  approximated  by  those  in  L  °  F  This  is  equivalent  to  the 

statement  that  for  any  T  €  ^  functions  in  $jlT  cannot  be  well 

approximated  by  those  in  M.  1|  ,  since  M, |  is  either  {0},  K  °  F"1, 

j— 1  1  j  1  R  T 

^  -1 

or  K  o  f  for  one  of  the  affine  transformations  F  :  T  T.  Applying 
this  observation  to  all  triangles  T  6  ^  yields  the  following  lemma, 

which  shows  that  functions  in  ^  cannot  well  approximate  those  in  . 
This  lemma  was  shown  previously  by  Bank  and  Dupont. 

A 

Lemma  4.10.  If  v  6  ^j-i  anc*  w  e  M  the  strengthened  Cauchy  inequaltiy 

(4.5.5)  a(v,  w)  <  y  ||| v (J)  |||w||| 

holds  for  fixed  y  e  (0,  1)  independent  of  v,  w  and  j  _>  2. 

Proof .  The  proof  will  be  done  locally  over  the  triangles  of  T^_^. 

Let  T  €=  Tj-i  *3e  arbitrary.  Then  there  exists  an  affine  transformation 

F_:  T  -*■  T.  Letting  w£  fl,  be  arbitrary  we  can  define  w*  w  °  F„.  Then 
1  K  J  1 

A  A  _ 

w*  £  K  or  w*  €  K  even  when  T  ?  x.  ,  and  w*  =  0.  Also  we  have 
k  b  j-i 

tv*  =  v  o  for  v  €  c  L  where  the  inclusion  may  be  strict 

because  of  boundary  conditions  or  interelement  continuity  conditions. 
Thus , 

Lnf  |v  °  F  -  w*|  .  21  inf  |v*  -  w*  |  .  >_  c  |w*|  . 

vGM  H  (Tr)  v*eL  H  (Tr)  H  (Tr) 


inf  || v  o  F  -  w*||  >_  inf 

L  (Tr)  v*GL 


,*  -  w*  II  2  >  C  ||w*  It  2 

L  (Tr)  L^(Tr) 


where  the  inequalities  on  the  right  are  from  lemma  4.8.  By  lemma  4.9 


these  inequalities  can  be  converted  into 
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inf  ]v  -  w 
v€EM 

J-l 


hV) 


>  c  w 


H1(T) 


-V! I|v  •  “lV> 


>  c  w 


l  2 

L  (T) 


where  the  constants  here  are  independent  of  d^,.  Combining  these 


relations  yields 


inf 


vein 


||v  -  w I)2 

IT-*- 


j-l 


_>  inf  |v  -  w|2  +  inf  |)v  -  w|f 

H'l  (T)  v^Af^_i  H  (T)  v€Mj_1  L  (T) 


>  c(|w|2  +  IMI2  )  =  c  II w || 

H  (T)  L  (T)  H  (T) 


Then  summing  over  the  triangles  of  Tj_^ 


inf  1 1  v  —  w  ||  1  =  inf  Z 

v£M.  ,  H  v€|U.  ,  T£T. 


j-l 


j-l  j-l 


II v  -  w If2 1 

H  (T) 


Z  inf  || v  -  wj 

Te-r.  .  veM. 

J-l  j-l 


I  i 

H  (T) 


>  c  Z  ||w II2  =  c  ||w ||2  . 

Ter.  .  h!(T)  h 

J-l 

So  using  the  equivalence  of  the  H'''  and  energy  norms 

inf  HI v  -  w|||  >  c  |||w|||  . 

^Vi 

This  is  equivalent  to  (4.5.5).  To  see  this,  assume  |||w|||  =  1.  Then 
inf  j 1 1 v  -  w|||  >  inf  |||v  -  w|||  >  c  |||w|||  =  c  , 

IIMII-i  veM 

v^M .  ,  -1 

J-l 


so  if  v 


w  =  1 


a(v,  w)  -  1  -  |  HI  v  -  w|||2  <  1  -  c  =  (1  -  c)  HI  v  HI  !||  w  HI  . 


The  general  case  of  (4.5.5)  follows  by  linearity.  □ 
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One  more  lemma  is  needed,  which  will  be,  for  our  purposes,  the 

A 

most  important.  It  states  that  functions  in  ,  the  kernel  of  the 

interpolation  mapping  from  to  Mj  are  all  quite  oscillatory. 

Intuitively,  the  relatively  smooth  functions  in  M,  agree  quite  closely 

with  their  interpolate  in  Only  the  most  oscillatory  functions  in 

M  can  have  the  zero  function  as  their  interpolate.  This  observation 

can  be  made  precise  quite  easily  through  the  use  of  lemma  4.9.  For  our 

purposes  a  function  is  "oscillatory"  when  its  norm  or  seminorm  is 

2 

much  larger  than  its  L  norm.  That  functions  in  are  oscillatory  in 
this  sense  is  shown  by  the  following  lemma. 


Lemma  4.11.  If  w  £  ii ,  ,  j  >  2  then 

(4.5.6)  7  II w[|  2  <  hj  |w|  x  <  c  || w||  2  , 

L  H  L 

with  c  >  1  independent  of  w  and  j. 

Proof.  It  is  enough  to  show  this  locally  on  the  triangles  of  x^  To 
see  this,  suppose  that  for  all  w  €  and  T  £  x^  ^ 

(4.5.7) 

with  c  >  1  independent  of  w,  T,  and  j.  Then  for  any  w  G  M 


7  IMI  2  1  h,  M  i  1  c  IMI  2 

L  (T)  J  H  (T)  L  (T)  ’ 


j 


2  IMI2 2  1  h4  Z  Mi  1  c2  Z 

c  T^X  L  (T)  3  T€T  H  (T)  T^X 

which  is  the  same  as  (4.5.6). 


ii  2 
L  (T) 


To  show  (4.5.7),  let  w  e  and  let  T  be  an  arbitrary  triangle 

in  x  .  Then  there  is  an  affine  transformation  F:  T  -*•  T  such  that 
j-i  R 

A  A 

w*  =  w  o  F  is  in  K  or  K  .  Since 
K  O 


,  is  a  norm  on  K  and  on  £  and 
.1  R  u 


f 


u  €  ^  or  u  €  K 
R  6 


7 IMI 2  i  lul  i  i c  IMI 2 

L  <Tr)  H  (Tr)  L  (Tr) 


with  c  >  1  independent  of  u.  Then 


7  ||w*||  <  |w*|  .  <  c  ||w*| 

^  *  *•  /  rf'  \  . »  -1-  /m  \ 


L  (Tr)  H  (Tr) 


1  2 

L2(Tr) 


so  using  lemma  4.9 


7  ||w||  2  <  h  |w|  <  c  ||w|| 

L  (T)  J  H  (T)  L  (T) 


completing  the  proof.  □ 


4.6.  Approximation 

In  this  section  the  approximation  inequality  (4.3.4)  needed 
to  justify  the  convergence  and  complexity  results  of  sections  4.3  and  4.4 
will  be  proven.  This  inequality  generalizes  a  result  of  Bank  and  Dupont 
for  quasi-uniform  grids.  Assume  the  elliptic  equation  is  H  regular 
for  a  £  (0,  1).  Then  their  result  is  that  if  is  a  family  of 

quasi-uniform  spaces  and  and  C7 ,  j  2,  are  defined  as  before,  the 
following  inequality  holds  for  r)  6  , 


III n  -  nil!  f  c  ej"  |||n|||  , 

where  r)  is  the  elliptic  projection  of  n  on  A  domain  with  a  slit, 

the  worst  case  ordinarily  arising,  yields  a  =  y  and  their  result 
becomes 


1 


’  "4 


9 


III n  -  SHI  <  C  e4/4  lllnlll  , 

the  same  as  the  result  here. 

The  result  here  generalizes  theirs  in  the  sense  that  it  applies 
to  locally  refined  grids  as  well  as  to  quasi-uniform  grids,  although  the 
exponent  of  0^  is  usually  poorer  in  our  inequality.  This  is  of  secondary 
importance  for  multi-level  convergence  since  0(N)  complexity  can  be  shown 
for  any  positive  exponent  of  0  ^ .  Only  the  hidden  constant  in  the  0(N) 
bound  is  affected  by  the  exponent. 

Enough  has  been  said  about  this  approximation  result  that  no 
further  introduction  seems  necessary.  The  technique  of  proof  is  quite 
interesting  and  begins  with  a  sequence  of  strengthened  Cauchy  inequalities 
that  follow  more  or  less  directly  from  the  results  of  the  last  section. 
Each  of  these  inequalities  is,  in  effect,  a  lower  bound  on  the  dihedral 
angle  between  two  subspaces. 


Lemma  4,12.  Let  w  G  ,  n  €  ,  s  €  and  o  e  0^  Then  the 

following  strengthened  Cauchy  inequalities  hold: 


(4.6.1a) 


a(n,  w)  <  c  0j  lllnlll  HI w|||  , 


(4.6.1b) 


a  (s ,  w)  <  c  0j^  (||  s  | 


(4.6.1c) 


a(n,  o)  <  C  0j/2  e~^2  III n | 


Proof .  For  (4.6.1a)  let  n  be  expanded  as 

J(e.) 

j  (i) 

n  -  2  ai  . 

n-1 
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Then  the  corresponding  finite  element  data  is 


J(6.) 

f  =  Z  a.  xfj)  , 

1  -|  '  1  * 


i=l 


11  *1 


with 

(4.6.2) 

Then 


a(ri,  <Ji)  =  (f,  4>)  ,  le  M.  . 


J(V 


2  t\ 0)\2 


i=l 


Thus 

(4.6.3) 


J(6.) 

i0jANd)  Z  “i  Alj)  H^j)H2 

J  j  i=l  1  1 

J(9.) 

=  9  ,  X<j)  ^  a*  ||H,<j)|||2 

J  j  i=l  1 


-  ei  \ 


fill  (9,  ^))1/2 
J  j 


Setting  ({)  =  w  in  (4.6.2) 

(4.6.4)  a(n,  w)  =  (f,  w)  <  |ff  ||  ||wf|  <  (0.  X<j))1/2 

J  j 


so  using  lemma  4.11, 
(4.6.5) 


■<n,  »)  <  C  h  (9  x')))1/2 

-  j  3  NJ 


<  c  8j 


using  the  eigenvalue  bound  (2.3.2).  The  proof  of  (4.6.1b)  is  the  same 
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« 


except  0  and  the  eigenfunctions  and  eigenvalues  {i|/j-l)}  j~l 
d  i  i=  1 

(j”l)  ^1“1 

^i=l  occur  throughout.  Thus,  instead  of  (4.6.5)  we  get 

*<■• »)  i  =  VV!  ^'j>)1/2  Ml  fll*lll 

1  «  p-1  Ills'll  IIHII 

HI  s|||  IIHH  , 

since  by  assumption  p  =  2.  For  (4.6.1c)  we  have  as  in  (4.6.4) 

a(o»  o)  <  (9j  xj  )1/2 

Expanding  o  in  eigenfunctions 


j-i 


Thus , 


j-1 

Z 

i-J(0j  1)+l 


llkp 


-1)  III  2 


j-1 


V2  .(j-D  lli,i(j-l)|l  2 


i-J(6j_1)+l 

i  i 

N, 

>  e  ,  x(j-1J 

J- 

E 

"  J'1  Nj-i 

1-J<eJ 

-  0  x(j_1) 

j-1  nj-i 

lloll2  . 

a2  ||^j' 


Dll  2 


and 


so 


» 
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a(n 


.  o)  £  (Sj  ^’)1/2  (Sj.!  ^j;"> 


-1/2 


By  the  eigenvalue  bound  (2.3.2) 

i  '  hi2  • 

j 

The  inequality  in  the  eigenvalue  bound  (2.3.2)  also  goes  the  other 

way.  One  way  to  see  this  is  to  expand  w  in  eigenfunctions  in  lemma  4.11: 

N. 

3 

w  =  I  a.  ib. 
i-1  1  1 

Then  the  left  side  of  (4.5.6)  is 

Nj  Nj  Nj 

7  II  z  '/'ill  1  h4  I  £  <*j  'M  1  c  h.  III  S  a  *  ||| 

i=l  3  i=l  1  „1  J  i=l  1 

H 

where  the  last  inequality  is  from  the  equivalence  of  the  H1  and  energy 
norms.  It  follows  easily  from  this  that 

^  ± ■ 

so  we  get 

a(n,  o)  <  C  0  e}/2  e-;(2  |[M||  |||o||| 

<  C  ej/2  ej2'2  III n III  lllolll  , 

since  the  mesh  ratio,  p,  is  fixed.  This  completes  the  proof.  □ 

Each  of  equations  (4.6.1)  of  this  lemma  states  that  the  energy 
inner  product  of  a  relatively  smooth  function  with  a  relatively  oscillatory 
function  is  small  compared  to  the  product  of  their  norms.  In  other  words, 
the  angle  between  these  functions,  as  measured  by  the  energy  inner 


^4 
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product.  Is  large.  In  the  appropriate  limit  as  9^  on  6^  ^  goes  to  zero 
each  of  equations  (4.6.1)  expresses  an  orthogonality  relation. 

Our  aim  is  to  show  that  the  function  in  can  be  well 
approximated  in  the  coarser  finite  element  space  ^  when  0^  is  small. 
Let  ^  be  the  orthogonal  complement  of  Mj  ^  in  with  respect  to  the 
energy  inner  product.  We  approach  the  question  of  approximation  of 
functions  of  in  M.  ^  by  the  back  door.  Instead  of  showing  directly 
that  M  ^  contains  a  function  near  any  function  in  S  ,  we  show  that 

3_i  3 

Sj  and  Mj  ^  are  nearly  orthogonal.  Since  consists  of  relatively 
smooth  functions,  by  analogy  with  the  lemma  just  proven  we  should  try 
to  show  consists  of  oscillatory  functions.  This  is  the  content 

of  the  following  lemma. 


Lemma  4.13.  If  u  £  ^j-1’  t^ien 

u  =  s  +  o  +  w  , 

A 

for  s  €  o  e  0^  and  w  G  and 


(4.6.7a) 


(4.6.7b) 

(4.6.7c) 


Mil  <cej!;  |||u| 

lolll  1  C  lllulll  , 
Ivlll  <  c  lllulll  , 


with  c  >  0  independent  of  u  and  j . 

Proof.  Since  ^  Q  ,  we  can  set 


u  ■  v  +  w  , 


for  some  v  €  w  6  Then  using  the  elementary  inequality 


xy  <  |  (x2  +  y2)  ,  x  ,  y  £  R  , 


and  lemma  4.5.3 


« 
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IIWII2  -  l||»  +  will2  -  jllvfll 2  +  IIHII2  +  2a(v,  w> 

>  IIMII2  +  IIWII2  -  2,  iiwii  IIMH 
1  a  -Y  HllWII2  +  lllwlll2)  . 


Thus , 

(4.6.8a)  III v|||  <  c  HHII  , 

(4.6.8b)  lllwlll  <  c  |||u|||,f 

where  c  >  0  depends  only  on  the  constant  y  €  (0,  1)  from  lemma  4.10, 
giving  (4.6.7c).  Now  since  u  =  v  +  w'€  ^  , 

a(v  +  w,  <(>)  =  0  ,  <t>  G  M  ^  , 


so  v  is  the  elliptic  projection  of  -w  on 


Writing 


v  =  s  +  o  , 

for  S  e  Sj_v  O  e  , 

a(s  +  w,  <t>)  *  0  ,  (J)  e  , 

a(o  +  w,  <J>)  »  0  ,  0  e  0  j  , 


since  ^  and  (7^  ^  are  energy  orthogonal.  That  is,  s  and  o  are  the 
elliptic  projections  of  -w  on  0 ^  respectively.  Then 


lllslll2  *  -a(s,  w)  <  c  ej^  |||s|||  HI w|||  , 
by  lemma  4.12,  (4.6.1b).  Thus,  using  (4.6.8b)  we  have 

lllslll  5  c  e222  IIWII  <  C  e}22  IIWII 

and 

illo  III  <  lllwlll  <  c  IIWII  , 


completing  the  proof.  □ 


« 


The  lemma  just  proved  shows  that  the  orthogonal 

complement  of  ^  in  ,  is  contained  in 

V  ©  °i-i  ©  "j 

l 

and  for  any  u  £  M.  the  component  of  u  in  S.  ^  is  relatively  small 

1 

when  ^  is  small.  In  other  words,  ^  consists  of  relatively 

oscillatory  functions.  Since  5^  consists  of  smooth  functions  when 

0.  is  small,  one  might  expect  that  the  inner  product  of  functions 

in  S.  with  functions  in  M.  ,  would  be  small  when  0.  was  small.  This 
J  j-1  j 

is  the  content  of  the  next  lemma. 


Lemma  4.14.  If  n  e  5^  and  u  ©  ^j_i>  then 

a(n,  U>  £  C  ejM  III n III  |||u|||  . 

Proof.  Since  0^  is  the  only  member  of  the  family  of  parameters 

occuring  either  implicitly  or  explicitly  in  the  statement  of  this 
lemma,  we  are  free  to  choose  the  remaining  parameters  arbitrarily. 

In  particular,  0^  ^  ©  (0,  1)  may  be  considered  a  free  parameter  until 
its  value  is  chosen  at  the  end  of  the  proof.  Let  r)  and  u  satisfy  the 
hypotheses.  According  to  lemma  4.13,  u  can  be  decomposed  as 


u  “  s  +  o  +  w  , 


with  s  €  S 


|_1>  0  e  0j_!  and  w  G 


.  Then 


a(n,  u)  -  a(r|,  s)  +  a (n*  °)  +  a(n,  w) 


Bounding  each  of  these  terms  separately 


a(n,  S)  <  IIMH  lllslll  <  c 


n  u 


by  lennna  4.13,  (4.6.7a).  By  lemma  4.12,  (4.6.1c), 


9 


art,  o>  i  c  ej22  e-_122  IHnlli  ilHH  <  c  ej22  |||n||[  |||„|||  , 
where  the  last  inequality  is  from  lemma  4.13.  Similarly  by  lemma  4.12, 


(4. 6. Id), 


art,  V)£C  e222  Hindi  IIMII  £C  e222  Hlnlll  ||l»lll  • 


Combining  these  inequalities 


a(n,  u)  <  c(ej^  +  eJ/2  0-i/2  +  eJ/2)  |||n|||  |||u| 


Choosing  0^  ^  =  0  completes  the  proof.  □ 


Our  first  approximation  result  is  a  simple  consequence  of 


this  lemma. 


Theorem  4.15.  Let  n  €  5^  and  let  r|  6  ^  be  the  elliptic  projection 

of  n  on  Mj 


i(n  -  n,  4>)  =  o  ,  <j>  e  M  . 


In  -  nil!  <  c  oj'"  |||n|||  . 


Proof.  We  have 


111 n  -  n I)) 2  -  a(n,  n  -  n)  5  c  ej/4  |||n|||  ||| n  -  n|||  , 

by  lemma  4.14,  since  ri  -  ri  £  Mj  Dividing  by  |  | r)  -  n|  |  the 
result  follows.  □ 


An  important  feature  of  the  result  just  proven  is  that  the 


constant  c occurring  here  depends  entirely  on  local  properties  of  the 


t 
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finite  element  spaces  ^  and  .  In  fact,  c  depends  only  on  the 
constants  in  the  eigenvalue  bounds  and  on  those  of  the  last  section 
on  interpolation.  As  a  consequence,  it  is  readily  computable. 

Moreover,  the  approximation  property  just  shown  does  not  depend 
significantly  on  either  the  domain  or  the  boundary  conditions.  This 
is  in  contrast  to  the  result  of  Bank  and  Dupont  where  both  the 
exponent  of  0  and  the  constant  occurring  depend  on  the  regularity 
of  the  PDE  and  hence  on  the  domain. 

The  approximation  result  just  proven  is  almost  the  result 
we  needed  in  section  4.3.  Its  defect  is  that  it  is  expressed  in 

N, 

terms  of  the  subspace  S.  and  indirectly  the  eigenfunctions  , 

J.v  N  ini 

rather  than  in  terms  of  and  indirectly  Since  it  would 

be  difficult  to  analyze  the  simultaneous  displacement  iteration  of 

Nj 

section  4.3  in  terms  of  the  eigenfunctions  we  must  prove 

the  analog  (4.3.4)  of  theorem  4.15  which  has  replaced  by  5^ .  Exactly 
the  same  thing  occurred  in  chapter  three.  While  it  may  be  possible 
to  get  the  required  analog  directly  from  theorem  4.15,  it  is  easier  to 
repeat  the  steps  in  its  proof  with  slight  modifications.  Actually, 
this  could  have  been  done  in  the  first  place,  but  it  seemed  more 
natural  this  way.  First  we  show  the  analog  of  lemma  4.12.  Only 
need  be  replaced  by  S ^ .  There  is  no  necessity  to  replace  ^  and 

0,  1  since  these  subspaces  play  only  an  intermediate  role. 


Lemma  4.16.  Let  w  6  Hj,  n  e  Sj ,  .s  6  Sj  ^  and  o  6  0^  Then  the 
following  strengthened  Cauchy  inequalities  hold: 


« 
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(4.6.9a)  a(n,  w )  <  c  G^2  j||n|||  |||w|||  , 

(4.6.9b)  a(s,  w)  <  c  G*72  |||s|||  |||w|||  , 

(4.6.9c)  a(n,  o)  <  c  G*72  6./72  ||tn|||  ||M||  • 

Proof.  Equation  (4.6.8b)  is  the  same  as  (4.6.1b)  so  is  already  established. 
For  (4.6.8a)  expand  n  in  the  eigenfunctions  {ipp^}, 

J(0.) 

n  =  £  ct.  , 

i=i 


and  let  f  be  given  by 


J(6 . ) 

3 

f  =  I 
i=l 


Then 


(4.6.10)  a(n,  <p)  =  b  (f ,  #)  ,  <J)  e  M  . 

Now  instead  of  (4.6.3)  we  get  by  a  similar  argument 

llflt  <  «,  *JJ)>1/2  lllnlij  , 

3 

so  with  <p  =  w  in  (4.6.10), 

(4.6.11)  a(n,  w)  =  b(f,  w)  <  (0  A,Jj))1/2  |||n|||  IM^  . 

Then  as  in  (3.3.3),  since  the  eigenvalues  of  the  mass  matrix  lie  on 
an  interval  f^,  o]  for  o  >  1  independent  of  j, 

(4.6.12)  llwllhia172  ||w||<  ca1/2  h.  |||w|||  , 

where  we  used  lemma  4.11.  Combining  (4.6.11)  and  (4.6.12)  and  using 


the  fact  that  0  is  fixed 


♦ 
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*<n.  ')!•  ye,  x«>),/2  III n III  IIMII 
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by  the  eigenvalue  bound  (2.3.2).  For  (4.6.9c)  we  have  as  in  (4.6.11) 

a(n,  O)  <  (0J  A^))1/2  Hlnlll  ||o|^  <  o1/2  (6.  A<j))1/2  III n III  ||o||  . 

Then  as  in  (4.6.6) 

IMli  <e  xO-1),-^  IHolll  , 

J  J-l 

SO 

a(n.  o)  <  o1/Z(9  <9  X<J)  r1/2  Hlnlll  |||o III 

j  J  j-l 

i  c(9j  h-2)1/2  <9J_1  h£)-l/2  Hlnlll  lllolll  , 

using  the  eigenvalue  bounds  as  before.  Then  since  the  mesh  ratio  is 
fixed , 

a(n*  o)  <  c  ej/2  e~^2  |||n|||  |||o||| 

as  required.  □ 


Since  we  are  retaining  ^  and  0 ^  ^  and  only  replacing 
by  Sj ,  there  is  no  need  to  alter  lemma  4.13.  Then  following  the  same 
proof  as  for  lemma  4.14,  but  using  the  strengthened  Cauchy  inequalities 
of  lemma  4.16  rather  than  lemma  4.12,  we  get  the  analog  of  lemma  4.14. 
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Lemma  A. 17.  If  n  e  5  and  u  £  M.  .  then 
-  J-1 

a(n,  u)  £  c  6j/4  lllnlll  lli  U  III  . 

The  required  approximation  theorem  follows  easily  from  this 
lemma  just  as  the  analogous  theorem  4.15  followed  from  lemma  A.1A. 


Theorem  A.  18.  Let  r|  €  and  let  r)  e  ^  be  the  elliptic  projection 


of  n  on  M .  ^ , 


a(n  -  n,  4>)  =  0  ,  0  e  M 


Then 


<4. 6.13)  llln  -  nlll  <  C  6»/4  lllnlll  . 

With  this  approximation  result,  the  convergence  and  complexity 
results  of  sections  A. 3  and  A. A  are  now  completely  justified.  As  for 
theorem  A. 15,  the  constant  c  occurring  here  depends  only  on  local 
properties  of  the  finite  element  spaces.  For  particular  spaces  one 
could  easily  compute  this  constant  and  determine  a  rigorous  upper  bound 
on  the  rate  of  convergence  of  the  multi-level  iterations  of  section  A. 3. 
This  generalizes  the  heuristic  local  Fourier  mode  analysis  to  irregular 
finite  element  grids,  and  is,  moreover,  completely  rigorous.  The  hidden 
constants  in  the  0(N)  complexity  bounds  of  section  A. A  are  not  so 
easily  computed  since  they  depend  on  the  generally  unknown  constant  in 
the  error  estimate  (A.A.l).  However,  for  most  purposes  it  is  enough 
to  have  a  bound  on  the  spectral  radius  of  the  multi-level  iteration,  which 
this  analysis  gives  since  this  is  what  determines  an  iteration’s 


practicality. 


« 
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